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Preface 

In  the  Air  Force  Research  Laboratory’s  Propulsion  Sciences  and  Advanced  Concepts  Divi¬ 
sion,  teams  of  scientists  and  engineers  are  engaged  in  research  aimed  at  making  rockets  and 
aircraft  more  capable  and  more  efficient.  They  design  new  high-energy  molecules  to  im¬ 
prove  performance  of  propellants.  They  develop  and  investigate  new  materials,  searching 
for  ways  to  make  lighter-weight  components  that  function  at  higher  temperatures.  They 
attempt  to  make  vehicles  and  weapons  more  survivable  by  manipulating  their  exhaust 
chemistry  and  probing  that  of  adversaries’  vehicles  and  weapons. 

Some  of  this  research  is  done  analytically;  most  experimentally.  In  many  cases,  the 
substances  involved  in  the  experimental  investigations  are  extremely  flammable,  explosive 
or  toxic.  To  work  with  such  substances  safely  is  a  difficult  and  expensive,  but  necessary 
proposition.  If  it  were  possible  analytically  to  predict  the  outcomes  of  reactions  and  the 
conflgurations  of  molecules,  the  need  for  such  laborious  experiments  could  be  reduced  to 
simple  one-time  verifications  of  the  computational  predictions.  For  now  and  for  the  fore¬ 
seeable  future,  however,  analytical  approaches  are  vastly  slower,  more  difficult,  and  more 
expensive  than  experimentation.  Great  strides  are  needed  in  both  numerical  algorithms 
and  computer  hardware  before  computer  calculations  can  feasibly  replace  even  relatively 
simple  chemical  experiments  on  the  types  of  molecules  most  valuable  in  aerospace  propul¬ 
sion  applications.  This  project  originated  in  this  need,  which  is  only  one  example  among 
many  beneficial  potential  applications  of  quantum  modeling  of  chemical  systems. 

No  great  strides  have  resulted,  but  perhaps  a  step  was  taken  in  the  right  direction, 
thanks  to  the  generous  support  of  numerous  people  and  organizations.  This  project  was 
sponsored  by  the  High  Energy  Density  Matter  Star  Team  at  the  Air  Force  Research  Lab¬ 
oratory  Propulsion  Directorate.  My  thanks  to  the  team  s  leaders.  Dr.  Jeff  Sheehy  and 
Dr.  Mario  Fajardo,  as  well  as  their  redoubtable  branch  chief.  Dr.  Pat  Garrick,  for  the 
opportunity.  Our  mutual  leaders,  Dr.  Steve  Rodgers  and  Dr.  Phil  Kessel,  and  recently 
Mr.  Mike  Huggins,  have  also  been  indispensably  supportive.  I  was  also  honored  to  receive 
encouragement  and  counsel  regarding  the  mechanics  of  completing  a  doctorate  from  three 
of  the  most  senior  leaders  of  the  Propulsion  Directorate:  Col  John  Rogacki,  Col  Wesley 
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Cox,  and  Dr.  Robert  Corley.  I  hope  to  return  the  favor  some  day  by  passing  on  their  wise 
advice  to  other  struggling  young  students. 

Computers  were  crucial  to  this  project-the  more  and  the  bigger  the  better.  Key 
portions  of  the  work  were  facilitated  by  the  use  of  the  Scientific  Visualization  Laboratory 
at  the  Aeronautical  Systems  Center  High  Performance  Computing  Major  Shared  Resource 
Center.  Thanks  to  Dr.  Jerry  Boatz  for  helping  me  surmount  some  difficulty  in  using  this 
facility.  Thanks  also  to  Dean  Wadsworth,  Roy  Hilton,  and  Alan  Kawasaki  for  their  help 
getting  me  access  to  computers  at  the  Air  Force  Research  Laboratory. 

The  members  of  my  dissertation  committee,  Dr.  William  Bailey,  Dr.  Larry  Burggraf, 
Lt  Col  Jeffery  Little,  and  Dr.  Anthony  Palazzotto,  and  my  advisor,  Dr.  Dave  Weeks,  not 
only  provided  essential  guidance,  they  were  graciously  accomodating  regarding  the  fits  and 
starts  in  the  progress  of  my  research,  dictated  by  the  demands  of  my  job  at  the  Air  Force 
Research  Laboratory. 

And  there  were  those  who  served  by  being  there,  and  being  themselves-my  friends 
and  family.  I  could  not  begin  to  express  my  appreciation  for  the  countless  demonstrations 
of  their  support  they  have  given  me  these  past  six  years,  and  the  years  that  led  up  to  them. 
First  among  them,  a  veritable  geyser  of  wisdom,  faith,  encouragement,  and  inspiration,  is 
one  who  shuns  credit  and  recognition,  but  will  always  get  it  from  me  nonetheless. 

Michael  J.  MacLachlan 
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Transmission  coefficient  of  the  as5Tnmetric  square  well  for  several 
choices  of  propagation  time  step  At  for  the  Mpller  states  only.  The 
correlation-function  step  of  the  channel-packet  process  uses  Ate  = 
1.0  •  10~^  atomic  units  in  all  cases . 


The  error  in  the  transmission  coefficients  of  the  asymmetric  square 
well  transmission  function  as  computed  in  the  Schrodinger  picture. 
Error  is  computed  using  equation  (5.24),  with  the  analytic  square- 
well  potential  as  the  reference . 

Convergence  of  the  transmission  coefficient  for  the  asymmetric  square 
well  as  the  coordinate  spacing  Ax  is  decreased  in  the  Schrodinger 
picture.  The  accuracy  of  the  discrete  representation  of  the  square 
potential  is  also  improving  with  decreasing  Ax.  Error  is  computed 
using  equation  (5.24),  with  the  analytic  square- well  potential  as  the 
reference . 

Error  in  transmission  coefficients  for  asymmetric  trapezoidal  well  po- 
tentals  (dotted  curve),  compared  to  those  for  the  corresponding  asymmtet- 
ric  square  wells  on  coarser  numerical  grids  (solid  curve).  The  trape¬ 
zoidal  wells  have  no  grid  points  within  the  well  walls  for  Ax  =  0.01, 
one  for  Ax  =  0.02,  three  for  Ax  =  0.04,  and  so  on,  as  enumerated 
in  Table  5.4.  The  entire  channel-packet  calculation  is  done  in  the 
Schrodinger  picture.  Error  is  computed  using  equation  (5.24),  with 
the  analytic  square-well  potential  as  the  reference . 

Convergence  in  the  Schrodinger  picture  of  the  transmission  function 
for  a  trapezoidal  well  of  fixed  shape  as  coordinate  spacing  Ax  is  varied. 
The  reference  transmission  function  uses  Ax  =  0.0025  atomic  units. 

All  calculations  use  At  =  1.0  •  10“®  atomic  units . 

Convergence  with  decreasing  time  step  At  of  transmission  functions 
computed  in  the  Schrodinger  picture  using  Mpller  states  produced  in 
the  interaction  picture,  compared  to  transmission  functions  produced 
entirely  in  the  Schrodinger  picture.  Error  is  computed  using  equation 
(5.24),  with  the  analytic  square- well  potential  as  the  reference.  .  .  . 
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for  trapezoidal  wells  modeled  in  the  Schrbdinger  picture  (dashed  line) 
is  included  for  comparison.  Error  is  computed  using  equation  (5.24), 
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6.3.  The  reactant  channel  packet  propagated  backward  in  time  to  f  = 

—2000  atomic  units,  and  the  product  channel  packet  propagated  for¬ 
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6.4.  The  reactant  Mpller  state,  (solid  line),  is  the  result  of  propagat¬ 

ing  the  intermediate  reactant  state  forward  in  time  to  t  =  0,  where  the 
interaction  and  Schrodinger  pictures  are  again  identical.  The  dashed 
line  is  the  product  Mpller  state,  the  result  of  propagating  the 

intermediate  product  state  backward  in  time  to  t  =  0 .  6-7 

6.5.  The  absolute  value  squared  of  the  correlation  function,  computed  us¬ 
ing  equation  (2.26) .  6-8 

6.6.  The  transmission  coefficient,  computed  from  Mpller  states  generated 

using  the  interaction  picture  (dotted  line)  and  the  SchrOdinger  picture 
(solid  line) .  The  interaction-picture  S-matrix  elements  were  computed 
on  a  grid  half  the  size  required  for  the  Schrodinger-picture  matrix 
elements .  6-9 

6.7.  The  fom  triple  Gaussian  potentials  used  to  investigate  the  effect  of 

varying  potential  slope  on  the  accuracy  of  propagation  in  the  interac¬ 
tion  picture .  6-11 

6.8.  The  growth  in  amplitude  error  (equation  (4.2))  during  the  computa¬ 

tion  of  the  reactant  Mpller  state  for  Potential  1.  The  solid  line  corre¬ 
sponds  to  a  scond-order  nested  interaction-picture  Lanczos  propaga¬ 
tor,  the  alternating  dots  and  dashes  to  a  first-order  one.  The  dotted 
line  shows  the  error  in  the  computation  of  the  same  Mpller  state  using 
a  split-operator  propagator  in  the  SchrSdinger  picture.  All  propaga¬ 
tions  shown  used  a  time  step  At  =  1.0  atomic  unit.  The  amplitude 
error  for  all  three  propagations  is  measured  against  a  split-operator 
propagation  using  At'  =  0.1  atomic  units .  6-12 


Xll 


Page 


Figure 

6.9.  The  growth  in  amplitude  error  (equation  (4.2))  during  the  computa¬ 

tion  of  the  reactant  Mpller  state  for  Potential  2.  The  solid  line  corre¬ 
sponds  to  a  scond-order  nested  interaction-picture  Lanczos  propaga¬ 
tor,  the  alternating  dots  and  dashes  to  a  first-order  one.  The  dotted 
line  shows  the  error  in  the  computation  of  the  same  Mpller  state  using 
a  split-operator  propagator  in  the  Schrodinger  picture.  All  propaga¬ 
tions  shown  used  a  time  step  At  =  1.0  atomic  unit.  The  amplitude 
error  for  all  three  propagations  is  measured  against  a  split-operator 
propagation  using  At'  =  0.1  atomic  units . 

6.10.  The  growth  in  amplitude  error  (equation  (4.2))  during  the  computa¬ 

tion  of  the  reactant  M0ller  state  for  Potential  3.  The  solid  line  corre¬ 
sponds  to  a  scond-order  nested  interaction-picture  Lanczos  propaga¬ 
tor,  the  alternating  dots  and  dashes  to  a  first-order  one.  The  dotted 
line  shows  the  error  in  the  computation  of  the  same  Mpller  state  using 
a  split-operator  propagator  in  the  SchrSdinger  picture.  All  propaga¬ 
tions  shown  used  a  time  step  At  =  1.0  atomic  unit.  The  amplitude 
error  for  all  three  propagations  is  measured  against  a  split-operator 
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Abstract- 

The  interaction  picture  is  used  together  with  the  channel-packet  method  in  a  new 
time-dependent  approach  to  compute  reactive  scattering  matrix  elements.  The  channel- 
packet  method  enables  the  use  of  the  interaction  picture  for  computing  reactive  S-matrix 
elements  by  splitting  the  computational  effort  into  two  parts.  First,  asymptotic  reactant 
and  product  wavepackets  are  individually  propagated  into  the  interaction  region  of  the  po¬ 
tential  to  form  Mpller  states.  The  interaction  picture,  in  contrast  to  the  usual  Schrodinger 
picture  of  quantum  mechanics,  is  so  constructed  that  a  wavefunction  that  experiences  no 
change  in  potential  (that  is,  a  free-particle  wavefunction)  remains  always  fixed,  with  no 
translation  or  distortion.  In  the  Schrbdinger  picture,  free-particle  wavefunctions  translate 
and  spread  with  time.  By  removing  free-particle  spreading,  the  use  of  the  interaction  pic¬ 
ture  reduces  the  size  of  the  region  of  space  that  must  be  modeled  when  computing  the 
Mpller  states.  Since  the  asymptotic  wavepackets  are  propagated  in  time  independently  of 
each  other,  it  is  possible  to  choose  an  asymptotic  Hamiltonian  and  corresponding  inter¬ 
action  picture  that  is  well  suited  for  each  arrangement  channel.  By  using  two  different 
interaction  pictures,  one  for  the  reactant  arrangement  channel  and  one  for  the  product 
arrangement  channel,  it  is  possible  to  realize  savings  in  the  required  grid  size.  During  the 
second  part  of  the  channel-packet  computation,  the  reactant  and  product  wavepackets  ob¬ 
tained  from  the  first  part  of  the  calculation  are  further  propagated  using  the  Schrodinger 
picture.  The  time-dependent  correlation  between  the  evolving  wavepackets  is  calculated 
as  they  split  into  energetically  accessible  arrangement  channels  and  are  absorbed  using 
absorbing  boundary  conditions. 

The  use  of  the  interaction  picture  for  computing  S-matrix  elements  is  developed, 
validated,  and  illustrated  using  a  simple  one-dimensional  reactive  example  where  the  size 
of  the  grid  required  for  computing  the  Mqller  state  in  the  interaction  picture  is  reduced  by 
a  factor  of  two  when  compared  with  required  grid  size  in  the  Schrodinger  picture.  Larger 
reductions  in  grid  requirements  are  realized  when  the  wavepackets  remain  compact  while 
evolving  into  Mpller  states,  especially  when  reactant  or  product  momenta  are  high. 
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Application  of  the  Interaction  Picture  to  Reactive  Scattering  in  One  Dimension 


/.  Introduction 

1.1  Context  of  the  Research 

The  purpose  of  this  project  is  to  advance  certain  computational  approaches  to  molec¬ 
ular  scattering  theory,  part  of  the  physics  that  underlies  chemical  kinetics.  Chemical  reac¬ 
tions  are  macroscopic  manifestations  of  numerous  individual  collisions  among  molecules  of 
various  species,  carrying  some  thermal  distribution  of  translational,  rotational,  vibrational, 
and  electronic  energies.  Each  pair  of  colliding  molecules  experiences  a  spatially  dependent 
mutual  force,  which  is  appreciable  over  some  “interaction  region”  of  relatively  small  inter- 
molecular  separation  distances  and  falls  to  zero  in  the  asymptotic  limit  of  large  separations. 
The  interaction  force  is  the  observable  manifestation  of  an  interaction  potential-energy 
function  whose  effect  is  significant  in  the  interaction  region,  and  which  approaches  a  con¬ 
stant  value  in  the  asymptotic  limit.  A  classical  view  of  a  molecular  collision  event  is 
illustrated  in  Figure  1.1.  The  figure  shows  an  elastic  collision  in  a  coordinate  system  in 
which  one  of  the  particles  is  fixed.  In  a  center-of-mass  coordinate  system,  comparable 
dynamics  could  occur  if  one  of  the  particles  were  of  much  larger  mass  than  the  other.  The 
light  particle  travels  in  a  straight  line  through  the  region  of  space  where  the  intermolecular 
potential  is  negligible,  along  the  asymptotic  paths  noted  in  the  figure.  The  interaction 
region  is  small  compared  to  the  asymptotic  regions  of  the  interaction. 

The  example  in  Figure  1.1  is  an  elastic  collision,  wherein  the  colliding  particles  are 
changed  in  momentum  only;  not  in  internal  configuration.  Over  the  course  of  more  com¬ 
plex  “reactive”  collisions,  munerous  types  of  changes  may  occur  in  the  colliding  reactant 
molecules,  including  the  rearrangement  of  atoms  within  or  between  molecules,  intercon¬ 
version  of  kinetic  and  potential  energy,  and  the  exchange  of  kinetic  energy  among  its 
translational,  rotational,  vibrational,  and  electronic  manifestations.  However,  the  event 
may  always  be  seen  as  consisting  of  a  brief  period  of  strong  interaction,  preceded  and 
followed  by  long  periods  of  negligible  interaction. 
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Figure  1.1  Schematic  of  a  classical  elastic  collision  event,  as  it  might  occur  between  a 
free  and  a  fixed  particle  in  a  laboratory  reference  frame,  or  between  a  very 
light  and  a  very  heavy  particle  in  a  center-of-mass  frame(l). 

In  order  to  predict  the  outcome  of  a  reaction,  or  to  understand  a  known  reaction 
from  first  principles,  one  would  like  to  know  the  probabilities  that  various  asymptotic 
post-collision  states  will  result  from  the  asymptotic  pre-collision  states.  The  informa¬ 
tion  needed  to  derive  the  probability  that  reactants  in  each  possible  energetic  state  and 
physical  configuration  will  give  rise  to  each  possible  state  and  configuration  of  products 
is  contained  in  the  scattering  operator  S.  The  state-to-state  reaction  probability  can  be 
averaged,  weighted  by  the  statistical  momentum  distribution  of  the  reactants,  to  obtain 
scattering  cross  sections,  which  can  be  used  in  turn  to  calculate  reaction  rates  and  pre¬ 
dict  macroscopic  reaction  outcomes.  The  scattering  operator  is  derived  numerically  in  an 
energy  representation  of  finite  order,  known  as  the  S  matrix.  The  square  of  the  absolute 
value  of  the  S  matrix  gives,  for  each  pair  of  reactant  and  product  energies,  the  probability 
that  a  reaction  yielding  the  particular  set  of  products  at  the  selected  product  energy  and 
configuration  will  result  from  the  collision  of  the  reactants  at  the  selected  reactant  energy 
and  configuration. 

Both  time-dependent  or  time-independent  techniques  have  been  developed  to  com¬ 
pute  the  S  matrix  (or  portions  of  it)  numerically.  Time-independent  methods(2-6)  are 
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limited  in  speed  and  scope  of  application  by  their  dependence  on  the  diagonalization  of 
large  matrices.  Time-dependent  methods(7-14)  have  been  hampered  by  their  requirement 
to  follow  the  motion  of  quantum  or  semiclassical  particles  over  large  changes  in  position 
and  momentum  during  the  course  of  the  collision,  over  the  full  time  from  before  the  re¬ 
actants  begin  to  interact,  through  the  collision,  to  the  time  when  the  products  no  longer 
interact  with  each  other.  Both  the  time-independent  and  the  time-dependent  methods 
therefore  present  large  computational  problems,  which  can  only  be  solved  with  reason¬ 
able  speed  and  accuracy  for  collisions  involving  very  few  atoms.  Based  on  the  numerical 
techniques  available  in  1988,  Kosloff  estimated  that  fast  minicomputers  could  reasonably 
be  used  for  problems  with  three  coordinate  degrees  of  freedom,  and  that,  based  on  ex¬ 
pected  improvements  in  computational  hardware  alone,  this  capability  would  increase  by 
one  degree  of  freedom  every  seven  to  eight  years(7).  The  analytical  techniques  Kosloff 
envisioned,  like  the  ones  treated  in  this  project,  all  invoke  the  Born-Oppenheimer  approx¬ 
imation.  This  is  the  assumption  that  electrons  adapt  instantaneously  to  the  motion  of 
their  nuclei,  allowing  the  treatment  of  atoms  without  dealing  with  electrons  as  separate 
objects(15).  Within  this  approximation,  the  neglect  of  overall  system  rotations  enables 
three  degrees  of  freedom  to  suffice  for  examining  two-atom  collisions,  as  will  be  seen  in 
Section  1.3.  Each  additional  atom  adds  three  more  degrees  of  freedom  to  the  problem, 
and  the  computational  time  required  to  compute  the  S  matrix  increases  rapidly  with  the 
number  of  degrees  of  freedom.  The  scaling  of  computational  effort  with  degrees  of  freedom 
is  estimated  in  Appendix  A.  Despite  this,  recent  algorithmic  developments  such  as  absorb¬ 
ing  boundary  conditions  have  already  yielded  some  six-dimensional  calculations  within  the 
eight- year  period  following  Kosloff ’s  estimate(  16-19).  This  is  a  good  illustration  of  the 
extra  leverage  algorithmic  advances  offer  to  complement  progress  in  computer  hardware. 
However,  computational  improvements  notwithstanding,  existing  computational  capabili¬ 
ties  are  still  far  from  sufficient  for  full  quantum-mechanical  treatment  of  the  majority  of 
reactive  scattering  problems. 

Typically,  quantum  computations  are  performed  using  a  view  or  “picture”  of  the 
mechanics  of  quantum  wavefunctions  known  as  the  Schrbdinger  picture.  Complex  wave- 
functions  containing  both  the  positional  and  momentum-space  distributions  of  a  quantum 
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system’s  probability  density,  are  common  to  all  quantum  “pictures.”  The  absolute  value 
squared,  (’I'|  'J'),  of  a  wavefunction  ’4’  gives  the  distribution  of  the  probability  of  the  sys¬ 
tem  being  found  in  a  particular  configuration.  The  Schrodinger  picture  is  characterized  by 
wavefunctions  obeying  the  Schrodinger  equation  of  motion  in  the  form 

=  (1.1) 

where  H5  is  the  full  system  Hamiltonian  as  discussed  in  Section  1.3.  Scattering  calculations 
can  be  performed  fruitfully  in  the  Schrodinger  picture,  but  this  project  demonstrates  that 
another  framework,  the  interaction  picture,  offers  computational  advantages  at  least  in 
some  scattering  problems.  In  the  interaction  picture,  described  in  detail  in  Section  3.5, 
wavefunctions  obey  the  Schrodinger  equation  in  the  same  form, 

1^)/  =  H/  1^)/  ,  (1-2) 

as  in  the  Schrodinger  picture,  but  the  Hamiltonian  operator  H/  and  the  wavefunction  ipj 
are  both  defined  in  such  a  way  as  to  cancel  out  all  “free”  motion  in  the  system — such  as 
the  motion  experienced  in  the  scattering  system  in  its  asymptotic  region. 

In  this  project,  the  interaction  picture  is  part  of  a  collection  of  techniques  assembled 
to  increase  the  efficiency  with  which  S-matrix  elements  can  be  computed.  Three  new  time- 
dependent  computational  approaches — the  channel-packet  method,  absorbing  boundary 
conditions,  and  the  interaction  picture — are  used  together  in  a  single  model  for  the  first 
time.  The  channel-packet  method  allows  the  calculation  of  relevant  sets  of  reactive  S- 
matrix  elements  without  the  expense  of  approximating  the  entire  matrix,  and  is  the  key 
to  surmounting  long-standing  difficulties  with  the  application  of  the  interaction  picture  to 
reactive  scattering  calculations  (1).  The  other  two  techniques  complement  the  channel- 
packet  method  by  decreasing  its  computational  effort.  The  assembly  of  techniques  is 
developed  theoretically  in  two  dimensions  in  Section  1.3,  and  analyzed  computationally 
in  one  dimension  elsewhere  in  this  docmnent.  It  is  extensible  in  a  straightforward  way  to 
systems  with  higher  degrees  of  freedom,  and  should  improve  further  in  speed  if  transported 
to  parallel-architecture  computers. 


1-4 


1.2  Research  Approach  and  Scope 

The  tools  developed  for  this  project  are  based  on  simulated  time  evolution  of  wave- 
functions  as  they  interact  with  one-dimensional  potentials.  Both  the  interaction  picture 
(Section  3.5)  and  the  Schrodinger  picture  are  used.  The  channel-packet  method,  described 
in  Chapter  II,  requires  two  separate  time  evolutions,  one  of  which  (to  calculate  Mpller 
states)  benefits  from  the  use  of  the  interaction  picture,  and  one  of  which  (the  calcula¬ 
tion  of  correlation  functions)  is  done  in  the  Schrbdinger  picture  with  the  aid  of  absorbing 
boundary  conditions  (Section  3.4).  The  split-operator  technique  of  Section  3.1  is  the  al¬ 
gorithm  of  choice  for  propagation  in  the  Schrodinger  picture,  which  is  used  as  a  reference 
for  comparison  with  the  interaction-picture  results.  Two  potentially  valuable  approaches 
to  wavepacket  propagation  in  the  interaction  picture  were  set  aside  when  they  failed  to 
produce  accurate  results  reliably.  These  methods,  the  finite-basis  approach  and  Tannor  s 
“Heisenberg”  approach,  are  described  in  Chapter  IV,  along  with  a  third  approach  that  is 
accurate  and  straightforward  to  implement,  but  cannot  provide  a  means  of  reducing  compu¬ 
tational  grid  sizes.  The  only  interaction-picture  technique  that  consistently  gave  accurate 
results  on  a  variety  of  potentials  is  the  sequential  nested  interaction  picture,  described  in 
Section  3.5.4.  Results  based  on  the  application  of  the  sequential  nested  interaction  picture 
appear  in  Chapters  V  and  VI. 

1.3  Quantum  Reactive  Scattering  in  One  and  Two  Dimensions 

The  research  reported  here  deals  primarily  with  one-dimensional  quantum-mechanical 
models  of  reactive  scattering.  To  illustrate  how  the  techniques  can  be  applied  to  problems 
of  higher  order,  the  two-dimensional  theory  is  also  developed.  Both  theories  are  developed 
in  the  familiar  Schrodinger  picture  of  quantum  mechanics. 

In  time-dependent  calculations,  scattering  events  are  modeled  using  complex  wave- 
functions  that  evolve  in  a  vector  space  that  is  viewed  in  two  equivalent  representations. 
In  the  position  or  “coordinate”  representation,  the  amplitude  of  the  wavefunction  gives 
the  distribution  of  its  probability  amplitude  across  the  position  coordinate  x.  In  the  mo¬ 
mentum  representation,  obtained  by  Fourier  transform  from  the  position  representation, 
the  amplitude  of  the  wavefunction  represents  its  momentum,  p  =  ^k.  It  is  often  conve- 
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nient  to  express  momentum  in  terms  of  k  rather  than  p.  This  is  done  freely  henceforth 
in  this  document,  particularly  since  the  system  of  atomic  units  is  always  used,  in  which 
Planck’s  constant  h  —  1.  A  system  involving  n  particles  of  mass  mj,  then,  is  described 
by  a  single  wavefunction  I'i/’)  of  reduced  mass  //,  position  (x)  =  (^|x|^),  and  momentum 
(p)  =  (V’l  p  \il’)  in  a  position  and  a  momentum  space  of  dimension  3(n  —  2).  The  reduction 
in  dimension  from  3n  to  3(n— 2)  is  accomplished  by  fixing  the  center  of  mass  and  neglecting 
any  overall  rotation(20). 

Time-dependent  methods  are  based  on  numerical  solutions  of  the  time-dependent 
Schrodinger  equation, 

ihj^\^it))=H\,p{t)),  (1.3) 

which,  for  time-independent  Hamiltonian  operators  H,  has  the  formal  solution, 

(1-4) 

where  \ip  (to))  is  the  state  vector  at  time  t  =  to-  Therefore,  the  Hamiltonian  conveniently 
expressed  as 


H  =  Ha  +  V(x),  (1.5) 

must  be  derived.  Given  a  numerical  approximation  for  the  interaction  potential  V (x)  be¬ 
tween  the  react, ants,  all  that  remains  is  to  add  the  asymptotic  potential  Ha,  containing  the 
kinetic  energy  of  the  reactants’  translation  toward  each  other,  plus  any  internal  vibrational 
or  rotational  energy  the  reactants  carry.  The  asymptotic  potential,  Ha,  plays  an  impor¬ 
tant  role  in  the  construction  of  Mpller  operators  (Section  2.2)  and  the  interaction  picture 
(Section  3.5).  The  numerical  approximation  of  the  full  Hamiltonian  is  used  to  derive  the 
time-evolution  operator  U(t,io)  =  which  governs  the  behavior  of  the  system 

throughout  the  collision  process.  Time-dependent  numerical  methods  use  approximations 
of  U(t',  t)  over  successive  short  time  intervals  At  =  t'-t  to  evolve  the  system  wavefunction 
over  the  required  time  period. 
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1.3.1  The  One- Dimensional  Interaction  Hamiltonian.  For  purposes  of  illustra¬ 
tion,  the  simplest  collision  to  describe  is  that  of  two  structureless  particles  with  no  angular 
momentum.  In  this  case,  the  motion  of  the  wavefunction  relative  to  center-of-mass  coordi¬ 
nates  is  one-dimensional.  Using  bold  type  for  operator  variables,  the  Hamiltonian  for  this 
one-dimensional  system  can  be  written 

H  =  Ha  +  V(rc),  (1.6) 

where  the  potential  operator  V(3;)  =  U(x)  describes  the  interaction  of  the  two  particles, 
and  the  asymptotic  Hamiltonian, 


Ha  =  lim  H 

“  s->±oo 

2/i  ’ 


(1.7) 


is  obtained  in  the  asymptotic  limit  where  lim  V{x)  =  0.  Here  k  represents  the  one- 
dimensional  momentum  operator  conjugate  to  x.  Equation  (1.6)  is  not  only  a  convenient 
way  of  accounting  for  the  kinetic-  and  potential-energy  contributions  to  the  full  Hamil¬ 
tonian;  it  will  also  serve  as  a  motivator  and  guide  to  the  construction  of  the  interaction 
picture  in  Section  3.5. 


1.3.2  Two-Dimensional  Interaction  Hamiltonians.  If  one  of  two  colliding  par¬ 
ticles  is  a  diatomic  molecule  consisting  of  atoms  labeled  A  and  B,  and  the  other  particle 
is  an  atom  labeled  C,  the  resulting  three-body  system  requires  3(3  2)  =  3  coordinates 

to  locate  each  atom  in  the  rotationless  center-of-mass  reference  frame.  If  the  angle  given 
by  the  vertices  A,  B,  and  C  is  fixed — ^for  example,  if  the  atoms  are  constrained  to  move 
along  a  single  line— then  the  resulting  collinear  system  requires  only  two  coordinates  to 
locate  the  atomic  positions.  One  possible  choice  of  coordinates,  called  bond  coordinates, 
is  illustrated  in  Figure  1.2,  where  the  coordinate  x  is  the  distance  between  atoms  A  and 
B,  and  y  is  the  distance  between  atoms  B  and  C.  Bond  coordinates  have  the  advantage 
of  generality,  describing  all  possible  arrangements  equally  well.  The  Hamiltonian  in  bond 
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A 


B 


C 


Figure  1.2  Bond  coordinates  for  three  collinear  bodies. 


(1.8) 


(1.9) 


mBmc 
ms  +mc 


(1.10) 


are  the  reduced  masses  of  the  two  pairs  of  adjacent  atoms,  and  via-,  tub,  and  me  are  the 
masses  of  the  respective  atoms(21).  In  the  asymptotic  limit  of  large  x  (or  y),  the  atom  A 
(C)  no  longer  interacts  with  the  diatom  BC  (AB)  and  V  {x,y)  —y  V  {y)  (or  V  {x)).  The 
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disadvantage  of  bond  coordinates  is  manifest  in  these  limits,  where  the  Hamiltonian  in 
equation  (1.8)  is  not  separable  like  equation  (1.6)  because  of  the  kinetic  coupling  term 

Using  an  alternative  choice  of  coordinates,  called  Jacobi  coordinates  (Figure  1.3),  the 
asymptotic  Hamiltonian  becomes  separable,  at  the  cost  of  using  different  Jacobi  coordi¬ 
nates  for  different  arrangements  of  particles  in  the  asymptotic  regions  before  and  after  the 
collision.  In  general,  the  three  atoms  A,  B,  and  C  may  be  arranged  in  any  of  four  ways: 

•  all  three  atoms  may  be  relatively  close  together  and  bound  or  strongly  interacting 
(ABC), 

•  all  three  may  be  relatively  far  apart  and  not  interacting  (A  +  B  -h  C), 

•  A  and  B  may  form  a  diatom  while  C  is  relatively  distant  (AB  -h  C),  or 

•  B  and  C  may  form  a  diatom  while  A  is  relatively  distant  (A  +  BC). 
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These  four  cases  are  termed  “arrangement  channels”  in  the  context  of  molecular 
collision  theory.  However,  as  far  as  a  collision  event  is  concerned,  the  first  case  cannot  be 
an  asymptotic  state  of  the  reactants  or  products.  It  must  be  either  a  temporary  transition 
state  in  the  midst  of  the  collision,  or  a  bound  state.  A  bound  state  cannot  be  the  state  prior 
to  collision,  since  there  are  no  free  particles  to  collide  with  one  another.  Neither  can  it  be 
the  state  following  the  collision  in  a  conservative  system  which  began  in  an  unbound  state. 
For  computational  simplicity,  it  is  also  possible  to  exclude  the  three-body  case  A  -t-  B  -t- 
C  by  a-sfinmiTig  that  the  system  begins  in  one  of  the  two  remaining  arrangement  channels 
and  restricting  the  energy  of  the  collision  so  as  to  make  complete  dissociation  energetically 
inaccessible.  The  interaction  then  becomes  a  one  with  only  two  possible  outcomes,  the 
arrangements  AB  -I-  C  and  A  -I-  BC. 

Each  of  the  two  arrangements  has  a  corresponding  set  of  Jacobi  coordinates.  In 
Figure  1.3,  the  two  arrangements  are  labeled  with  the  subscripts  7  =  1  and  7  =  2.  In 
each  arrangement,  the  coordinate  Ry  represents  the  distance  separating  the  two  members 
of  the  diatom,  and  is  the  distance  between  the  free  atom  and  the  center  of  mass  of  the 
diatom.  The  Hamiltonian  in  Jacobi  coordinates  is  given  by 


(1.11) 


where  is  the  mass  of  the  lone  atom,  and  fj,^  is  the  reduced  mass  of  the  atom-diatom 
system  in  the  arrangement  channel  labeled  by  7.  Thus,  for  one  labeling  of  the  two  ar¬ 
rangements, 


= 


{rriA  +  tub)  me 
rriAmBmc 


(1.12) 


and 


f^2  = 


niA  (tub  +  me) 

mATUBTric 


(1.13) 
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In  the  asymptotic  limit  of  large  the  Hamiltonian  reduces  to  the  sum  of  a  free-particle 
Hamiltonian  and  an  internal  Hamiltonian, 


lim  HL 

^±oo 

2mj  2^1^ 

2m^  2//^ 


+  lim  V^(r^,iiy) 

r^— >>±00 

+  V7nt(i^) 


(1.14) 


In  equation  (1.14),  \]^^{Ry)  =  V.Y(r.y,ity)  is  the  internal  potential  of  the  diatom, 

and  kr.^  and  are  the  momentum  operators  conjugate  to  and  Ry  respectively.  The 
last  term  in  equation  (1.14), 

=  (1.15) 


is  the  Hamiltonian  of  an  isolated  diatom,  and 


H 


7  _ 

rel 


2m-y 


(1.16) 


describes  the  behavior  of  a  free  particle  of  mass  m.y  with  an  energy  equal  to  the  relative 
kinetic  energy  between  the  atom  and  diatom.  In  general,  ^  my,  Hy,  and 

V'^{ry,Ry)  ^  V'^{ry,Ry>)-,  the  two  Hamiltonians  are  as  distinct  as  the  arrangements 
themselves.  Thus,  in  an  exchange  reaction,  such  as  AB  +  C  — >  A  +  BC,  one  must 
use  both  coordinate  systems  in  order  to  take  advantage  of  the  simplified  form  (1.14)  for 
both  the  reactant  and  product  asymptotic  Hamiltonians.  The  requirement  for  multiple 
coordinate  systems  and  Hamiltonians,  and  the  need  to  change  coordinate  systems  in  mid¬ 
collision,  present  the  principal  difficulties  with  the  use  of  Jacobi  coordinates  in  scattering 
calculations. 


1-11 


II.  Time-Dependent  Moleevlar  Scattering  Theory 

The  channel-packet  method  is  a  technique  for  obtaining  selected  elements  of  the  S  matrix 
by  means  of  time-dependent  quantum  propagation  techniques.  Such  techniques  are  based 
on  the  numerical  solution  of  the  time-dependent  Schrbdinger  equation, 

ihj^\^l^it))  =  K\rP{t)).  (2.1) 

The  formal  solution  to  the  time-dependent  Schrbdinger  equation  for  a  time-independent 
Hamiltonian  operator  H  is 


l^Pit))  =  V{t,to)\^{to)) 

=  (2-2) 

where  {ip  (to))  is  the  state  vector  at  time  t  =  to,  and  U(t,io)  =  e-***(*-*'>)/^  is  the  time- 
evolution  operator.  Commonly,  it  is  convenient  to  set  to  =  0,  and  use  the  shorthand 
notation  U(t)  =  U(t,0). 

There  have  been  several  applications  of  the  channel-packet  method  to  a  variety  of 
scattering  problems.  In  particular,  the  method  has  been  used  to  compute  state  to  state 
S-matrix  elements  for  the  collinear  H  +  H2(n)  ^  H2(n')  +  H  reaction(22,23)  and  more 
recently  for  a  two  dimensional  model  OC  +  OH(n  =  0)  OCO(n  =  0)  +  H  reaction(24). 
In  a  larger  three-dimensional  calculation,  Dai  and  Zhang  have  used  the  channel-packet 
method  to  compute  exact  state-to-state  S-matrix  elements  for  the  H  +  O2  reaction(12). 
One  common  advantage  shared  by  all  of  these  calculations  is  the  facility  with  which  the 
channel-packet  method  provides  S-matrix  elements  over  a  wide  range  of  energies.  In  addi¬ 
tion  to  exact  quantum  calculations,  the  channel-packet  method  has  also  been  useful  in  for¬ 
mulating  several  approximate  strategies  for  computing  S-matrix  elements.  These  include 
a  new  semiclassical  method  for  computing  S-matrix  elements  developed  by  Garashchuk 
and  Tannor(25),  and  the  application  of  the  multiconfiguration  time-dependent  Hartree 
approach  to  computing  S-matrix  elements  by  Jackie  and  Meyer(26). 
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2.1  Channel  Packets 


The  two-dimensional  asymptotic  Hamiltonian  of  equation  (1.14),  H2,  has  two  compo¬ 
nents,  and  H/nj.  Since  both  and  have  known  eigenvalues  and  eigenstates, 
eigenstates  of  the  asymptotic  Hamiltonian  can  be  constructed  in  Jacobi  coordinates  as 
direct  products  of  free-particle  eigenstates,  =  I/C7),  and  diatomic  eigenstates,  I7): 


H2  1^7,7) 


2^  1^7)  It)  +  m)  It) 


Ifey’T), 


(2.3) 


where  the  index  7  is  now  used  to  label  the  channel,  a  term  that  includes  both  the  arrange¬ 
ment  channel  and  the  internal  eigenstate  of  the  diatom.  For  each  arrangement  chaimel, 
there  exists  one  channel  for  each  eigenstate  of  the  internal  Hamiltonian 

Each  channel  is  associated  with  a  single  eigenstate  I7)  of  the  diatom,  and  infinitely 
many  free-particle  eigenstates  1^7).  Rather  than  attempting  to  consider  the  non-square- 
integrable  states  |A:7)  I7)  =  1^7,7)  individually,  it  is  useful  in  collision  modeling  to  construct 
localized  linear  combinations  of  the  |A:7,7), 


^7n(o«t))  ~  J  ^  dk.yr]^  (k-y)  \k^,  7)  •  (2-4) 

Equation  (2.4)  defines  two  wavepackets,  or  “channel  packets,”  in  the  channel  labeled  7. 
The  channel  packet  is  expanded  in  terms  of  the  eigenstates  |fc7,7)  of  th®  7-chaimel 
Hamiltonian,  where  the  expansion  coefficients  7;^  (^7)  are  chosen  to  be  appreciable  over 
some  particular  range  of  momenta  ^7.  The  packet  |'02ut)  also  expanded  in  terms  of 
the  eigenstates  1^7,7),  with  expansion  coefficients  r]_  (k^).  The  channel  packets  and 
li’Zut)  represent  reactants  and  products,  respectively,  in  the  channel  7. 
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It  is  convenient  to  choose  the  expansion  coefficients  77^  {k~f)  so  that  the  state 
is  a  complex  Gaussian  in  the  momentum  coordinate  k-y, 


V±  (^7) 


_1 

27r  ( Afc)Q 

^  exp 

(^k'y  ^-yo) 


"h  (^7  ^70)  ? 


(2.5) 


where  ( Afe)Q  is  the  uncertainty  in  the  momentum,  and  the  constants  r^o  ^70  fix  cen¬ 
ter  of  the  Gaussian’s  representation  in  coordinate  and  momentum  space,  respectively(27). 
The  Gaussian  is  chosen  because  it  can  be  manipulated  analytically,  and  can  be  defined 
to  include  whatever  particular  range  of  values  of  momentum  may  be  of  interest.  Since 
the  Fourier  transform  of  a  Gaussian  is  also  a  Gaussian,  the  coordinate  representation  of 
such  a  channel  packet  in  Jacobi  coordinates  is  the  direct  product  of  a  Gaussian  in  the 
coordinate  and  a  diatomic  eigenfunction  in  the  Ry  coordinate.  Such  channel  packets  are 
the  basis  for  the  method  described  in  Section  2.3  for  evaluation  of  matrix  elements  of  the 
scattering  operator. 


2.2  The  Scattering  Operator 

The  scattering  operator  S  is  the  keystone  of  quantum  molecular  scattering  theory. 
The  scattering  operator  relates  the  reactant  and  product  states,  \ipin)  and  \'^out)^ 
interaction  as 


\4’ont)  =  S\4’in)-  (2-6) 

For  the  scattering  operator  to  exist,  the  interaction  potential  V (x)  must  satisfy  the  con- 
ditions(l) 

1.  Y(^)—0  as  a;  00;  that  is,  the  potential  approaches  zero  more  rapidly  than 

does  x~^  in  the  asymptotic  limit; 

2.  Y{x)=0  as  a;  — >  0;  that  is,  the  potential  becomes  unbounded  less  rapidly 

than  x~^  at  the  origin;  and 

3.  V(x)  is  continuous  nearly  everywhere  for  0  <  a:  <  oo. 
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It  is  important  to  note  that  even  though  the  pure  internuclear  Coulomb  potential  does 
not  meet  conditions  1  and  2  above,  the  presence  of  atomic  electrons  serves  to  make  the  S 
operator  formalism  applicable  to  atomic  and  molecular  problems.  Atomic  electrons  shield 
nuclear  charges  from  one  another  at  large  separations,  and  add  to  the  interatomic  repulsion 
at  small  distances. 

It  is  useful  to  express  the  scattering  operator  as  the  product 

S7'7=Q7;tf^7,  (2.7) 

where  the  M0ller  operators  are  defined  as 

al=  lim  (2.8) 

The  M0ller  operators  are  isometric  =  1),  and  are  unitary  only  if  no 

bound  states  exist.  Mpller  operators  also  obey  the  “intertwining”  relation(28) 

HfiX  =  (2.9) 

where  H  is  the  full  Hamiltonian,  and  H2  is  the  asymptotic  Hamiltonian  in  arrangement 

channel  7.  The  intertwining  relation  will  be  important  in  Chapter  3.  In  terms  of  the 
Mpller  operators,  the  probability  of  scattering  from  a  given  reactant  state  to  a  given 
product  state  given  by 

=  (2-10) 

where  7  and  Y  indicate  the  reactant  and  product  channels,  respectively.  The  states  It/’-  )  = 
and  1^/+)  =  IV'Zn)  are  called  Mpller  states. 

The  scattering  operator  is  customarily  evaluated  in  its  energy  or  momentum  repre¬ 
sentation,  called  the  S  matrix.  S-matrix  elements  are  expressed  in  the  momentum  repre- 
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sentation  as(31) 


(By  -  £,)  =  (<:;„7'|S|fcy7) 

=  Ii„7) 

=  {k'^',i-\k^,l+)  ■ 

The  density  of  states, 

dE 
“  d\k^\ 
d 

d\k.y\ 

n?\k^\ 


(2.11) 


(2.12) 


is  introduced  in  equation  (2.11)  to  convert  the  S  matrix  from  energy  normalization  to 
momentum  normalization.  The  absolute  value  squared  of  the  S-matrix  element  is 

the  probability  that  a  reaction  that  starts  with  reactants  in  the  state  1^7,7)  will  yield 
products  in  the  state  ky,  7'^ 


The  right-hand  side  of  equation  (2.11)  is  the  inner  product  of  the  vectors 


and 


1^:7, 7+)  —  1^7)7)  (2' 13) 


1^7' >7'-)  =  I  ^7' >7')  ,  (2.14) 
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which  are  eigenstates  of  the  full  Hamiltonian.  This  can  be  shown  using  the  intertwining 
relation  (2.9),  together  with  equations  (2.13)  and  (2.3): 

Hlfe.„7+)  = 

=  niKZ\k^,7) 

A  similar  argument  shows  that  the  states  \ky,'y-),  are  also  eigenstates  of  the  full  Hamil- 
tonian. 

Two  observations  regarding  the  states  |fc>y,7ib)  are  useful  in  the  development  of  the 
channel-packet  method  in  Section  2.3.  First,  1^7, 7±)  satisfy  the  orthogonality  relation 

()t;„y±|fc„7±>  = 

=  S^fy  (  kyt  ,  7  I  kyj  7  ^ 

=  SyyS  {k’y,  -  ky)  .  (2.16) 

Second,  if  the  full  Hamiltonian  H  is  time-independent,  then  solution  of  the  time-dependent 
Schrodinger  equation  (2.1)  for  the  time  evolution  of  |A:7,7±)  yields, 

\ky,'y±;t)  =  lJit)\ky,'y±) 

=  exp[-i(i|;fc?(,,,  +  E,M)|]|fc,.7±>.  (2.17) 

2.3  Computation  of  S-Matrix  Elements  using  Channel  Packets 

Channel  packets,  introduced  in  Section  2.1,  were  developed  and  used  by  Weeks  and 
Tannor  for  computing  matrix  elements  of  the  scattering  operator (29,30).  In  this  time- 
dependent  method,  a  pair  of  wavepackets,  |V’7„)  and  defined  using  equation  (2.4), 
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is  used  to  derive  S-matrix  elements  8^,’"''^^.  for  a  broad  range  of  momenta,  and  fcy ,  within 

7'’  ^ 

a  particular  pair  of  reactant  and  product  channels  7  and  7'. 

The  calculation  begins  by  applying  a  numerical  time-evolution  scheme  like  those 
developed  in  Sections  3.1  and  3.2  to  the  channel  packets  to  produce  the  M0ller  states 

(2.18) 

=  lim  \4>l„) 

t-^—oo 

«  ^-^-rlh^TXlrlh  1^7^^ 

and 

i,t)  = 

(2.19) 

where  ±r  are  values  of  time  t,  defined  such  that  essentially  all  of  the  product  and  reactant 
wavepackets  experience  only  the  asymptotic  potential  H„  for  all  t  greater  than  +t  and  less 
than  — r.  The  coordinate  representations  of  the  states  |V'7n)  chosen  to  be  in 

the  interaction  region  of  the  potential  at  time  t  =  0.  The  effect  of  the  M0ller  operator 
on  the  initial  state  |^7ra)  propagate  it  backward  from  time  t  =  0  to  time  t  =  — r  using 
the  asymptotic  Hamiltonian  H2  until  the  resulting  intermediate  state  exits  the  interaction 
region,  and  then  to  propagate  the  intermediate  state  forward  to  time  t  =  0  again  using 
the  full  Hamiltonian.  Similarly,  the  product  state  is  propagated  forward  from  time 

t  =  0  to  time  t  =  +t  under  H2',  and  back  to  time  t  =  0  under  H.  The  resulting  Mpller 
states  can  then  be  propagated  under  the  full  Hamiltonian,  using  the  methods  described  in 
Chapter  III. 

The  Fourier  transform  of  the  time-evolved  Mpller  state  (t))  =  |V’7n)  is 

/OO 

dt exp  (iBt)  U(t)  \iil) 

-OO 

roo 

=  dt  exp  {iEt)lJ {1)^1  lip] J  .  (2.20) 
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The  channel  packets  and  be  expanded  as  in  equation  (2.4),  in  terms  of  the 

eigenvectors  and  fey, 7'^  respectively.  Substitution  of  equations  (2.18)  and  (2.4) 

into  (2.20)  yields 


/oo  roo 

dt  exp  (iEt)  U(t)Q+  /  dk^t]^  (k^)  \k^,  7) 

-OO  J —OO 

/oo  roo 

dt  exp  (iEt)  U(t)  /  dk^t]^  (k^)  \k^,'y+)  ■ 

-OO  J—oo 


(2.21) 


Equation  (2.17),  substituted  in  equation  (2.21),  gives 


(E)')  =  J  dtexp^^^e  J  dkjT]^{ky)\kry,‘j+) 

=  /_”  ifexp  (iff)  exp  |^-|  (!=■,)  1*^.7+) 

=  27r  y  dk^6  {k^  -  A:^)^  v+  {k'r)  \kj,n/+) 

=  /  dk^  [6  (fey  -ky)+6  [ky  +  k^)]  ??+  {hy)  \ky,  7-I-) 

^  1^71  J-oo 

=  T^{??+(+fcy)|-fA^y,7-|-)4-T/+(-fey)|-fcT,,7-|-)},  (2.22) 


dtexp  ^ 


27r^^ 
h\k-y\  ^ 
2txix^ 
h  |A:..),| 


where  the  two  degenerate  eigenstates  corresponding  to  positive  and  negative  momenta  +ky 
and  —k-y  are  the  only  terms  of  the  integration  singled  out  by  the  delta  functions.  The  imier 
product  of  this  expression  with  the  product  Mqller  state  is  then 


i’l  AliE 


27r 

j  t^y‘ 

’A‘7 

h  \ 

1  \kj>\ 

\\k,\ 

+  7/1  {k 

27r 

/ 

-/^7 

T\ 

1  |A:y 

\ky\ 

(k’y)v^  (+fe,)<fc;„y-|  fc„7+) 
+  1)1  (fcy)')+(-<7)(*4'.7'-|*+.7+)  } 

gj|j{,i(+ow').(+Msri?.,7^ 

+  V-  {-ky)  r]+  i+k^) 

+  rf-  (+^7')  'n+  {-k'y)  S+fc7^ 

+  V-  (“^y)  V+  {-kj)  • 


(2.23) 
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The  channel-packet  method  allows  a  piecemeal  evaluation  of  the  S  matrix  by  consid¬ 
ering  four  special  cases  of  equation  (2.23),  where  the  expansion  coefficients  are  chosen 
so  that  either 


1.  7?+  (+^)  =  rf-  (~^7')  =  0,  with  (r.y  (-r))  >  0  and  (r'^,  (+t)^  >  0; 

2.  77+  i+k^)  =  rj*.  j  =  0,  with  {r-,  (-r))  >  0  and  (r'y  (-l-r)^  <  0; 

3.  77+  {-k-y)  =  77*  (-fcy)  =  0)  with  {r^  (-r))  <  0  and  (r'y  {+t)'^  >  0;  or 

4.  77+  (-k^)  =  77*  =  0,  with  {r^  (-r))  <  0  and  (r'^  {+t)'^  <  0. 

In  case  1,  for  example,  the  reactant  wavepacket  starts  to  the  right  of  the  origin,  with  all 
free-particle  momentum  directed  to  the  left.  The  product  wavepacket  finishes  up  to  the 
right  of  the  origin,  traveling  to  the  right.  Equation  (2.23)  reduces  to 


(2.24) 


in  this  case.  Referring  to  equation  (2.20),  it  can  be  seen  that  the  left-hand  side  of  equation 
(2.24)  may  also  be  written, 


dt  exp  {iEt)  U (t)  | 


V’’ 


+ 


=  J  dt  exp  (iEt)  (ipZ.  U(t)|'0Xy 


(2.25) 


This  expresses 


as  the  Fourier  transform  of  the  correlation  function. 


cp''y{t)=  (V’l'|u(i)|v'T). 


(2.26) 


The  correlation  function  (2.26)  is  the  inner  product  of  one  Mpller  state  and  the  time 
evolution  of  the  other  Mpller  state.  Both  Mpller  states  can  be  calculated  directly  by  the 
propagation  methods  of  Chapter  III,  as  applied  to  equations  (2.18)  and  (2.19).  Using 
equation  (2.26),  one  Mpller  state  is  evolved  in  time  from  t  =  —00  to  t  =  —00,  and  its  inner 
product  with  the  other  Mpller  state  is  evaluated  as  a  function  of  time.  Then  the  Fourier 
transform  is  performed  to  produce 


(E)y  By  rearranging  equation  (2.24),  one 
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quarter  of  the  S-matrix  for  the  energy  range  of  interest  is  then  given  by 


(2.27) 


The  other  three  cases,  2  to  4,  deliver  the  rest  of  the  S-matrix,  according  to  the  general 
formula, 


(2.28) 


providing  the  remaining  three  combinations  of  reactant  and  product  momenta. 
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III.  Computational  Methods 

The  channel-packet  method  requires  time-dependent  propagation  of  wavepackets,  first  to 
arrive  at  the  Mpller  states,  then  to  evaluate  the  correlation  function.  A  number  of  methods 
have  been  developed  (34)  to  solve  the  time-dependent  Schrodinger  equation  (2.1)  numer¬ 
ically  for  Two  such  methods,  the  split-operator  and  the  Lanczos  algorithms,  are 

discussed  in  Sections  3.1  and  3.2.  Additional  numerical  techniques  to  apply  these  methods 
more  efficiently  to  channel-packet  calculations  are  presented  in  Sections  3.4  and  3.5. 


3.1  Split- Operator  Propagation 

The  split-operator  scheme  of  Feit  and  Fleck  originated  in  a  homologous  problem 
in  optics(35,36).  The  technique,  suitable  for  use  with  time-independent  Hamiltonians, 
is  based  on  the  separation  of  the  time-evolution  operator  into  separately-diagonalizable 
kinetic-  and  potential-energy  portions: 


U(to  +  ^t) 


e-*(T+v)At//iu(fo) 

g-iVAi/fig-iTAi/ng-i[T,V]AiV2;i^U  (to)  +  0{At^) 


(3.1) 


where  T  and  V  are  kinetic-  and  potential-energy  operators  respectively,  and  [T,  V]  = 
TV  —  VT  is  their  commutator.  The  Zassenhaus  formula(37) 


gA+B  ^  gAgBg§[B,Al  ^  ^  [a,B]]) 


(3.2) 


is  required  in  order  to  obtain  equation  (3.1),  since  T  and  V,  in  general,  do  not  commute. 
The  split-operator  formulation  of  the  time-evolution  operator, 


U(to  +  At)  =  e-*VAt/2ftg-iTAt/fig-iVAf/2fiu  ^  q  ^^^3^  ^  (3  3) 

can  be  derived  from  a  Taylor-series  expansion  of  equation  (3.1).  Equation  (3.3)  has  the 
benefit  of  retaining  the  second-order  accuracy  of  (3.1)  while  involving  only  the  two  opera¬ 
tors  T  =  P^/2/i  and  V  by  themselves,  eliminating  the  commutator. 
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The  kinetic-energy  evolution  operator,  is  diagonal  in  the 

momentum  representation,  and  the  potential-energy  evolution  operator,  e~*VAt/2fi^  jg 
agonal  in  the  coordinate  representation.  Thus,  the  effect  of  the  time-evolution  operator 
on  a  state  \ijj  (to))  over  a  short  interval  At  can  be  calculated  easily  with  the  help  of  two 
fast  Fourier  transforms  (FFTs).  The  Fourier  transform  and  its  inverse  have  the  effect  of 
converting  vectors  back  and  forth  between  the  momentum  and  coordinate  representations. 
The  calculation 


Ip  (x,  to  -f-  At)  =  (xl  U(At)  (to)) 

=  (a;|  e-*VA</2fie-iTA</ftg-iVAt/2fi 

=  (xj  e-*VAt/2fi  f  l^./^  f  dk  \k)  (A:|  e-»TA«/fi  f  ^k' \k')  {k'\ 


X  J dx"  \x")  {x"\ e-»VAt/2fi  ! dx'"  \x'")  {x'"\  ^  (to)) 

J  dx' \x')  j  dk{x'  Ifc)  I  dk' 


=  J  dx'  {x\  e-*VAt/2fi  1^/^  J  1^^  I  g-iTAt/fi  1^/^ 
X  Jdx"{k'  \x")  j  dx'"  {x"\e-^^^^/^^\x'")7P{x'",to) 

=  J  dx'b  [x  -  x')  e-*'"(^)^*/2A  j  dk  {x'  \k)  I  dk'6  {k  -  k') 

X  j  dx"  {k  \x")  e-*'^(-")A*/2/i  J  dx'"6  {x"  -  x"')  v>  {x"',  to) 


-iV(x)At/2h_ 


V 

:  J  dx"( 


hi 


—ikx"—iV{x 


')^‘/2V(x",to) 


can  be  performed  using  an  iV-element  complex  vector  to  approximate  any  function  x  t)  = 

{x  |x(i))?  or  its  Fourier  transform. 


{k  IxW) 

-=  /  dxe~^'^'''x{x,to)  • 

V27r  J-oo 


(3.5) 


The  diagonal  operators  ^^d  ^an  be  represented  by  iV-element  arrays. 

At  each  time  step,  the  wave  function  is  multiplied,  element  by  element,  by  the  diagonal 
operator  elements,  e-iV{xi)At/2h ^  resulting  vector  transfornied  to  the  momentum 
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representation  by  FFT.  The  vector  is  then  multiplied  by  the  diagonal  operator 
transformed  back  to  the  coordinate  representation  by  inverse  FFT,  and  multiplied  by  the 
operator  por  a  time-independent  Hamiltonian,  the  time  evolution  of 

IV’  (t))  can  be  calculated  by  repeated  applications  of  equation  (3.4),  using  the  same  diagonal 
operators. 


3.2  Lanczos  Propagation 

The  split-operator  technique  is  a  fast,  robust  method  of  obtaining  solutions  to  the 
time-dependent  Schrodinger  equation  for  time-independent  Hamiltonians.  (34)  However, 
there  are  cases,  as  in  the  interaction  picture  (see  Section  3.5),  where  the  form  of  the  time- 
evolution  operator  makes  it  impractical  to  separate  the  evolution  operator  into  portions 
diagonalizable  by  FFTs.  In  such  cases,  the  wavepacket  may  be  propagated  by  the  Lanczos 
method(44,45).  In  the  Lanczos  approach,  the  Hamiltonian  is  represented  using  basis  vec¬ 
tors  of  a  Krylov  subspace  generated  by  the  Hamiltonian.  The  basis  vectors  of  the  Krylov 
subspace  are  formed  by  repeated  operation  of  the  Hamiltonian  on  some  initial  state  |V>q): 

\g,J  =  H-\<t>o).  (3.6) 


It  is  generally  neither  necessary  nor  desirable  to  use  more  than  a  few  basis  vectors.  Com¬ 
monly,  the  dimension  M  of  the  Krylov  subspace  is  between  5  and  8.  An  orthonormal  set  of 
basis  vectors  |g„)  is  obtained  from  the  non-orthogonal  |^„)  using  Gram-Schmidt  orthogo- 
nalization,  and  used  to  form  an  M  x  M  finite-basis  approximation  of  the  Hamiltonian(45), 


— >Hm  = 


/02  0  .  0 

(^2  a2  0  •  •  •  0 

0  /^s  as  '  ■  •  ‘  : 

:  0  ••  ••  Pm-1  0 

•  •  '•  Pm-1  Q!M-i  Pm 

0  0  •••  0  OM 


(3.7) 
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Matrix  elements  in  equation  (3.7)  are  given  by 

an  =  (gni  H  \qn) ,  (3.8) 

Pn  \qn)  =  H  |9n_i)  -  a„-i  |g„_i)  -  /3;_1  \qn-2) ,  (3.9) 

Pi  =  0,  ki)  =  ko) ,  ko)  =  I  } ,  (3.10) 

where  |  )  is  the  null  vector.  The  orthonormality  of  the  basis  vectors  \qn)  is  used  to  derive 
the  subdiagonal  values, 


=  ||H  kn-l)  -  an-l  kn-l)  -  Pn-l  kn-2)|l  ,  (3.11) 

from  equation  (3.9),  employing  the  norm  defined  by 

lllx)ll  =  ^/(x^x>.  (3.12) 

and  choosing  all  /3„  to  be  real,  since  employment  of  the  norm  renders  the  phase  arbitrary. 
Within  the  reduced  dimensionality  of  the  Krylov  subspace,  the  Hamiltonian  is  easily  di¬ 
agonalized.  Once  the  Hamiltonian  is  diagonalized,  the  time-evolution  operator  is  readily 
expressed  as  a  matrix,  and  the  time-evolved  state  vector, 

|V>(t  +  At))  =  e-®^*/'^|V’(t)),  (3.13) 

may  be  expanded  in  terms  of  the  normalized  basis  vectors  kn)-  This  process  is  iterated  to 
advance  the  initial  state  (to))  over  a  sequence  of  small  time  steps,  in  order  to  arrive  at 
a  final  state  (t)). 
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Suppose  the  propagation  is  performed  in  the  momentum  representation.  Then  the 
evolution  from  time  t  =  tj-i  to  t  =  tj  is  expressed, 

{k  {tj))  =  {k  |U  {tj,tj-i)  IV’  (tj-i)) 

M  M 

n=l  m=l 

Here  denotes  the  nth  Krylov  vector  in  the  basis  used  to  represent  the  Hamiltonian 
H(tj)  and  generate  the  state  vector  \'ip{tj))‘  The  choice  =  |V’(tj-i)))  using  the 
previous  time  step’s  state  vector  to  seed  the  present  time  step’s  basis,  allows  equation 
(3.14)  to  be  simplified  to 

n=l  m=l 

M  M 

n—\  m=l 

M 

n=l 

3.3  Computational  Bottlenecks 

The  channel-packet  method,  based  on  time-dependent  propagation,  provides  an  effi¬ 
cient  alternative  to  time-independent  close-coupled  type  calculations  (2-6)  for  computing 
S-matrix  elements.  However,  peculiarities  of  the  FFT  and  the  time-dependent  behavior  of 
Gaussian  wavepackets  both  may  cause  the  method  to  require  unnecessarily  large  coordinate 
grids  and  attendant  long  computation  times. 

The  FFT  is  much  faster  to  calculate  than  the  equivalent  discrete  Fourier  transform 
for  large  numbers  N  of  points  in  space  at  which  functions  are  evaluated.  However,  FFTs 
conflate  the  origin  xq  and  the  Nth.  grid  point  xn,  artificially  imposing  periodicity  on  the 
potential.  A  wavepacket  which  attempts  to  propagate  off  the  edge  of  the  grid  at  x  =  x^ 
instead  “wraps  around”  and  appears  on  the  other  side,  encoimtering  the  potential  at 
a;  =  0  instead  of  the  correct  potential  at  quasi-infinite  x.  The  same  problem  occurs  in  the 
momentum  representation  for  wavepackets,  the  absolute  value  of  whose  momentum  begins 


3-5 


to  exceed  the  value  that  defines  the  positive  or  negative  edge  of  the  momentum  grid. 
The  “wrapping  around”  of  all  or  part  of  the  momentum  representation  of  the  system’s 
wavefunction  results  in  immediate  severe  degradation  of  the  accuracy  of  the  simulation  as 
the  momentum  abruptly  reverses  and  takes  on  the  opposite  sign.  The  common,  brute-force 
way  to  avoid  these  grid-related  errors  is  to  make  the  grids  so  large  that  the  wavepackets 
never  approach  the  edges  in  either  representation.  The  additional  overhead  requirement 
is  compounded  by  the  fact  that  the  fastest  FFTs  constrain  the  allowed  number  of  grid 
points  N.  (In  the  simplest  implementation,  the  constraint  is  iV  =  2",  where  n  is  a  natural 
munber.) 

Another  factor  tending  to  force  the  grid  to  be  larger  is  the  well-known  spatial  spread¬ 
ing  of  Gaussian  wavepackets  governed  by  free-particle  Hamiltonians  (27).  For  the  channel- 
packet  method,  this  means  that  grids  must  be  large  enough  to  accommodate  spreading,  as 
well  as  translation,  of  the  product  and  reactant  states  as  they  are  propagated  backward 
and  forward  in  time  under  the  asymptotic  Hamiltonian.  Large  grids  lead  in  turn  to  long 
computation  times  for  the  FFT.  Both  of  these  grid-enlarging  factors  can  be  ameliorated, 
and  the  grid  size  reduced  significantly,  through  the  use  of  absorbing  boundary  conditions 
and  the  interaction  pictiue. 


3.4  Absorbing  Boundary  Conditions 

Absorbing  boimdary  conditions  can  combat  the  spurious  periodicity  of  the  FFT  grid 
by  including  an  imaginary  component  in  the  potential  V (x)  for  values  of  x  near  the  edges 
of  the  grid(23,38-43).  A  complex  potential  V{x)  =  V  (x)  ±f/(x),  for  some  real  function 
/  (x)  that  is  zero  over  most  of  the  grid  and  begins  to  increase  near  the  edges  of  the  grid, 
adds  a  real  exponential  decay  to  the  forward  and  reverse  time-evolution  operators.  Careful 
choice  of  the  function  /  (x)  can  assure  that  the  wavepacket  is  absorbed  before  it  crosses 
the  edge  of  the  grid,  and  is  not  reflected  or  transmitted  to  any  significant  degree. 


The  success  of  the  application  of  absorbing  boundary  conditions  to  the  channel- 
packet  method  lies  in  the  computation  of  the  correlation  function  (2.26).  The  product 
Mpller  state  ^  is  localized  in  the  interaction  region  near  the  origin,  so  the  correlation 
function  (t)  =  (ipZ.  U(t)j  1/’+^  need  be  evaluated  only  in  this  restricted  region,  once 
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the  M0ller  state  has  been  calculated.  Computing  the  correlation  on  an  FFT  grid  confined 
to  this  region  saves  time,  but  leads  to  invalid  results  if  absorbing  boundary  conditions  are 
not  employed  to  keep  portions  of  the  evolving  reactant  state  (t))  which  propagate  out 
of  the  interaction  region  from  “wrapping  around”  and  spuriously  returning  from  the  other 
side.  If  the  imaginary  potential  is  confined  to  the  portion  of  the  grid  where  the  product 
Mpller  state  is  zero,  it  does  not  effect  the  correlation  function,  so  S-matrix  elements  may  be 
derived  correctly  despite  the  altered  potential.  Calfas  and  Weeks(23)  have  demonstrated 
this  technique  in  two  dimensions  with  absorbing  boundary  conditions  of  the  form  /  (x)  — 
A  exp  (x  -  xo)^  /B^  for  the  collinear  H  +  H2  ^  H2  +  H  reaction,  using  the  Liu-Siegbahn- 
Truhlar-Horowitz  (LSTH)  potential (62-64). 

3.5  The  Interaction  Picture 

Propagating  the  wavepackets  in  the  interaction  picture  can  nearly  eliminate  spread¬ 
ing,  as  has  been  confirmed  for  simple  one-dimensional  potentials(46-49).  The  interaction 
picture  has  also  been  applied  to  molecular  predissociation  in  two  and  three  dimensions, 
where  a  single  arrangement  channel  and  corresponding  Jacobi  coordinates  are  sufficient  (50- 
52).  The  channel-packet  approach  allows  the  interaction  picture  to  be  extended  for  the 
first  time  to  reactive  scattering  requiring  two  arrangement  channels. 

3.5.1  Essentials  of  the  Interaction  Picture.  In  the  usual  (Schrodinger)  picture 
of  quantum  mechanics,  the  time-dependent  behavior  of  state  vectors  is  governed  by  the 
time  evolution  operator  U(t,to)  =  e"***^*"*®^/^,  for  time-independent  Hamiltonians  H. 
An  alternative  approach,  known  as  the  interaction  picture  (also  called  the  intermediate 
or  Dirac  picture),  originated  in  time-dependent  perturbation  theory.  In  time-dependent 
perturbation  theory,  the  Hamiltonian  operator  is  defined  as  the  sum  H  =  Hq  -I-  Hi,  where 
Ho  is  a  time-independent  Hamiltonian  susceptible  of  direct  analysis,  and  Hi  is  a  small, 
time-dependent  perturbation(53).  The  construction  of  the  interaction  picture  begins  with 
such  a  splitting  of  the  Hamiltonian(54). 
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In  the  interaction  picture,  the  Hamiltonian  is  most  commonly  expressed, 

H  =  Ho+V,  (3.16) 

where  Ho  is  the  asymptotic  Hamiltonian  and  V  is  the  potential,  as  in  equation  (1.5). 
However,  this  is  not  the  only  possible  expression  for  the  Hamiltonian,  and  is  not  necessarily 
the  most  useful.  Therefore,  let  us  retain  the  more  general  notation, 

H  =  Ho  +  Hi,  (3.17) 

where  Hq  is  any  time-independent  portion  of  the  full  Hamiltonian.  If  the  time-evolution 
operator  corresponding  to  Hq  is  denoted, 

Uo  (t)  =  (3.18) 


the  transformation  from  the  SchrOdinger  picture  to  the  interaction  picture  is 

=  (3.19) 

where  the  subscript  I  labels  the  interaction  picture,  and  the  subscript  S  refers  to  the 
Schrodinger  picture.  It  is  seen  readily  from  equation  (3.19)  that  lip  (0))j  =  ['^(O))^. 

The  Schrodinger  equation  in  the  interaction  picture  is 

=  -U+  (t)  Ho  \iP  {t))s  +  it)  (Ho  +  Hi)  1^  {t))s 
=  uJ(t)HiUo(t)|V’(i)>7 

=  (3.20) 


1^  (*))/ 
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where  the  interaction-picture  Hamiltonian  is  defined  as 


H/(t)-uS(t)HiUo(t).  (3.21) 

In  other  words,  H/  (t)  is  constructed  in  such  a  way  as  to  make  the  form  of  the  SchrSdinger 
equation  in  the  interaction  picture  the  same  as  in  the  Schrodinger  picture.  Hence,  the 
formal  solutions  of  the  Schrodinger  equation  in  the  interaction  picture  take  the  analogous 
form. 


=  (3.22) 


for  some  time  evolution  operator  U/(t,to)  initial  state  |•0(^o))/.  However,  since  the 
Hamiltonian  H/(t)  is  time-dependent,  the  simple  expression  U/(t,to)  = 
is  not  valid.  Instead,  integration  of  the  Schrodinger  equation  with  the  initial  condition 
U/  (^O)^o)  =  1  gives  the  perturbation  expansion. 


U/(t,to)  = 


1  +  ^  j*  Hi  (t')  Ui  (t'  -  to)  dt' 

1  +  ^  r  H/  (t')  dt'  +(^y  f  t  Hi  {t')  Hi  (t")  dfdt'  -\-  ■  ■  -(3.23) 

J to  \  ^  /  Jto  Jto 


so  called  because  its  most  common  application  lies  in  time-dependent  perturbation  theory. 
In  perturbation  theory,  the  time-dependent  portion  of  the  Hamiltonian  is,  by  construction, 
small  enough  that  the  series  (3.23)  converges  within  the  first  few  terms.  For  the  interaction 
Hamiltonian,  the  truncated  perturbation  series  may  not  be  of  sufficient  accuracy  for  a  given 
time  interval  At  =  t-  to-  The  Magnus  expansion, 

U,  (t,  („)  =  !+  (:^)  f  H,  ((')  dt'  +  i  (^)'  £  £  [H,  ((') .  H,  ((")]  dt"dt'  +  . . .  , 

(3.24) 


is  a  way  to  redress  the  convergence  problem(55,56).  The  commutator. 


[Hi  {t') ,  Hi  [t")]  =  Hi  (t')  Hi  {t")  -  Hi  (t")  Hj  (f)  (3.25) 
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in  equation  (3.24),  can  be  eliminated  by  using  the  equivalent  expansion(57) 

U,(t,(o)  =  {(')<#'  (3-26) 

-,4  f  f  (H'  (*')  H/  (*")  -  2e  («"  -  f)  H,  ((")  H;  {t'))dt’df  +  ..., 

Jto  .ho 

where  0  (t"  —  t')  represents  a  Heaviside  step  function.  In  time-dependent  numerical  meth¬ 
ods  employing  small-enough  time  steps  At,  first-order  truncations  of  the  series  (3.23), 
(3.24),  and  (3.26)  are  equivalent  and  accurate.  If  the  numerical  propagator  employs  time 
steps  that  are  large  enough  to  require  a  second-order  truncation,  equations  (3.24)  and 
(3.26)  provide  the  better  approximation. 

The  relationship  between  the  interaction-  and  Schrodinger-picture  time-evolution  op¬ 
erators  follows  from  the  definition  of  the  interaction-picture  state  vector  (equation  (3.19)), 

=  uj  (t)  U  (t,  to)  Uo  (to)  ItA  (to)) /  •  (3.27) 

Hence,  the  time-evolution  operator  in  the  interaction  picture  can  be  written 

Vi  (t,  to)  =  Uj  (t)  U  (t,  to)  Uo  (to)  •  (3.28) 

3.5. 1.1  Alternative  Formulations  of  the  Interaction  Picture.  The  construc¬ 
tion  of  the  terms  Ho  and  Hj  is  worth  a  moment  of  consideration.  Since  the  time  evolution 
of  a  state  in  the  Schrodinger  picture  occurs  as 

{tp  {t))g  =  exp  [-i  (Ho  -f  Hi)  (t  -  to)/h]  \ip  (to)) ,  (3.29) 

much  of  the  spreading  effect  of  the  asymptotic  Hamiltonian  is  counteracted  in  the  time 
evolution  of  the  interaction-picture  state  if  the  choice  Ho  =  Hq  is  made.  This  is  why  the 
interaction  picture  is  usually  constructed  this  way.  Another  attractive  choice  is  Ho  =  T  = 
/2fii) ,  the  kinetic-energy  operator,  which  simplifies  the  calculation  of  the  interac- 
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tion  Hamiltonian  H/  =  by  making  Ho  and  Hi  fimctions  of  momentum 

and  position  alone.  The  same  advantage  accrues  if  Hq  =  Hrei  =  h^k^/2/x.y  is  chosen.  The 
three  choices,  Hq  =  Ha,  Ho  =  T,  and  Ho  =  Hj-ei,  are  equivalent  in  one  dimension,  where 
Ha  =  T  =  Hrei. 

Two  additional  considerations  affect  the  choice  of  Ho.  First,  if,  as  in  perturbation 
theory,  the  effect  of  Hi  is  small  compared  to  that  of  Ho,  larger  time  steps  may  be  taken 
in  numerical  propagations  in  the  interaction  picture(49,50).  This  favors  Ho  =  Ha,  putting 
as  much  of  the  Hamiltonian  into  Ho  as  possible.  Second,  however,  the  interaction-picture 
methods  that  are  the  most  successful  in  reducing  computational  grid  size  (the  nested  in¬ 
teraction  pictures  described  in  Section  3.5.4),  are  implemented  using  Lanczos  propagation, 
which  depends  on  diagonality  of  Ho  in  the  momentum  representation  in  the  construction 
of  Krylov  basis  vectors.  In  more  than  one  dimension,  therefore,  the  best  choice  of  terms  is 
generally  not  Ho  =  Ha.  This  issue  is  discussed  further  in  Section  3.5.3. 

3.5.2  Scattering  in  the  Interaction  Picture.  The  Mpller  states  illustrate  the 
application  of  the  interaction  picture  to  time-dependent  scattering.  Since  Mpller  states  are 
defined  in  the  Schrodinger  picture  at  time  t  =  0,  they  are  equal  to  their  counterparts  in 
the  interaction  picture: 


\i’±)i  =  \^±)s- 


(3.30) 


For  the  same  reason,  the  asymptotic  reactant  and  product  states  are  likewise  invariant: 


Therefore,  for  example. 


^  in{out)  ^  j 


^in{out)  ^  g  ' 


(3.31) 


1^+)/  ~  |^-l-)s  —  ^+sl^7n)s 

=  (t,  0)  U2  (t) 

=  U^[t(0)  lim  V{0,t)VZ{t)\i>l)j 
«  Ut(0)U(0,-T)U2(-r)|^7J,,  (3.32) 


3-11 


where  the  unit  operator  Uq  (0)  is  introduced  to  suggest  the  form  of  equation  (3.28).  The 
general  form  of  the  Mpller  operators  in  the  interaction  picture  is 


nlj  =  ^  Um^  uS  (0)  U  (0,  t)  U2  (t) .  (3.33) 

Here,  our  choice  of  the  operator  Hq  comes  into  play.  If  the  selection  Hq  =  H2  is 
made,  then  U2  (t)  =  Uq  (t) ,  and  equation  (3.33)  gives 

QIj  =  ^nm^Uj(0)U(0,t)Uo(t) 

=  lim  U/(0,t).  (3.34) 

Another  potentially  useful  alternative  is  to  choose  Hq  =  sense  of  equation 

(1.16).  Since  and  commute,  equation  (3.33)  can  be  written, 

nlj  =  ^Um^U$(0)U(0,t)Uo(t)U7„,(t) 

=  Jm^UK0,t)U7„,(t),  (3.35) 

where  (t)  =  exp  The  action  of  these  Mpller  operators  on  channel  packets 

|^7n(emt))  Mpller  states 

\^±)l  =  ^±l\^ln{out)) 

= , 

«  .  (3.36) 

Thus,  with  the  appropriate  choice  of  Ho,  the  Mqller  operators  become  simple  propagations 
in  the  interaction  picture,  possibly  with  a  phase  shift. 

3.5.3  Wavepacket  Propagation  in  the  Interaction  Picture.  Numerical  propagation 
in  the  interaction  picture  is  not  a  trivial  problem,  considering  the  form  (3.21)  of  the 
Hamiltonian.  Interaction-picture  state  vectors  do  not  lend  themselves  to  split-operator 
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propagation.  Consider  a  simple  first-order  split-operator  type  propagation  over  a  single 
time  step,  choosing  Hq  = 


IV’ (*0 


+  At))j  ^  + 


IV’  ih)) 


=  exp » {(o)) 


=  exp 


f  iAt 

1 


,i/i(to+At/2)k^/2n 


xe 


-tfi(to+Ai/2)k2/2m. 


I IV’ (io)) 


(3.37) 


=  exD  -  I  _  !^gift(to+At/2)k?/2m^ 

^  \  2fi^  h 

H,)  e-i^(‘o+At/2)k2/2m^  1 1^ 

xV^  {r^,  R^)  e-*'^(‘o+At/2)k?/2m^  yiAtm^l/An^  1^  _ 


The  outer  operators,  are  directly  diagonalizable,  but  the  inner  exponential 

is  not.  Alternative  formulations  of  the  interaction  picture  fail  to  eliminate  this  problem. 
Thus,  beneficial  use  of  the  interaction  pictme  requires  a  finite-basis  approach  like  the 
Lanczos  method  of  Section  3.2,  whereby  the  evolution  operator  U/  (r)  may  be  diagonalized 
as  a  single  entity. 

The  Lanczos  approach  is  workable,  though  cumbersome  unless  further  refinements 
are  added.  The  computation  of  the  Krylov  basis  vectors  via  equation  (3.9)  involves 
evaluating  the  vectors 


H/  (t,)  14)/  =  uS  (t,)  HiUo  (tj)  14)/ ,  (3-38) 


which,  if  taken  head-on,  implies  propagation  of  the  state  vector  under  Hq.  The  effect  of 
equation  (3.38)  is  to  convert  4^^  to  the  Schrodinger  picture,  operate  upon  it  with  Hi, 
and  convert  the  result  back  to  the  interaction  picture.  If  Hq  is  not  a  free-particle  type 
Hamiltonian,  each  conversion  back  and  forth  involves  a  propagation  over  the  full  time  tj, 
not  just  the  incremental  time  At.  What  is  worse,  this  double  propagation  must  be  done  M 
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times  for  each  time  step  in  order  to  obtain  all  of  the  Krylov  vectors.  Therefore,  formulations 
of  the  interaction  picture  wherein  Hq  is  not  a  function  of  momentum  alone  incur  far  too 
much  computational  overhead  to  be  useful,  even  when  compared  to  Schrodinger-picture 
propagations  on  much  larger  grids.  For  free-particle  type  Hq,  propagation  is  not  necessary 
in  the  formation  of  the  basis  vectors,  but  the  implicit  conversion  to  the  Schrodinger  picture 
means  no  reduction  in  grid  size  is  achieved,  since  the  full  spreading  effect  of  Uq  (t)  is 
ultimately  incurred  before  being  reversed. 

The  direct  “sequential”  approach  just  described  to  computation  of  Krylov  vectors 
makes  Lanczos  propagation  in  the  interaction  picture  much  slower  than  propagation  in  the 
Schrodinger  picture  using  the  same  grid  and  time  step.  One  possible  way  around  this  ob¬ 
stacle  is  the  finite-basis  approach,  which  spreads  wavepackets  only  over  short  time  intervals 
and  does  not  use  Fourier  transforms.  The  finite-basis  approach  is  a  new  idea,  developed 
early  in  the  course  of  this  research  project,  which  failed  to  find  a  stable  computational 
implementation.  The  theory  and  some  computational  results  of  the  finite-basis  approach 
are  described  in  Section  4.3.  A  better-developed  idea  that  has,  in  contrast,  proved  com¬ 
putationally  robust,  is  the  nested  interaction  picture. 

3,5,^  Nested  Interaction  Pictures.  Tannor  and  others  have  observed  that  a  re¬ 
duction  in  the  number  of  required  grid  points  in  the  interaction  picture  can  be  achieved 
by  optimizing  the  grid  width  in  momentum  space,  or  the  grid-point  spacing  in  coordi¬ 
nate  space(49).  They  named  their  invention  the  nested  interaction  picture  because  the 
Hamiltonian  becomes  sandwiched  between  one  or  two  additional  pairs  of  operators. 

Two  forms  of  nested  interaction  pictures  have  been  developed:  the  so-called  P- 
adapted  and  PR-adapted  varieties.  In  the  P-adapted  picture,  the  state  vector  is 

1^)'/  = 

^  ^  (3  39) 
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where  the  expectation  value  of  the  momentum  operator, 


(P)  =  (V’l/P/ltA)/ 

=  |V')5 

=  (V>|5Ps|V’)s,  (3-40) 

is  taken  to  be  constant  for  the  duration  of  a  single  time  step.  The  unitary  transformation 
gi(P)Rs/ft  jjgg  effect  of  shifting  the  state  vector  in  such  a  way  as  to  make  its  average 
momentum  zero.  The  P-adapted  state  vector  obeys  the  wave  equation 

=H'|V»)^,  (3.41) 

where 


H'  =  (R5)  (3.42) 

The  PR-adapted  picture  shifts  the  origins  of  both  the  momentum  and  the  coordinate 
representations,  defining  the  state  vector, 

IV’)/  = 

^  gi(R>Ps/fig-i{P>Rs/fteiHo</ft|^^^^  (3_43) 

where  (R)  =  (V'l/  R/  |i/’)/  =  (V'b  ^5  \i’)s  expectation  value  of  the  position  operator. 
The  PR-adapted  equation  of  motion  is  given  by 

=  (3.44) 

where 


JJ//  _  gt<R)Ps/ftg-»(P>Rs/ftgiHot/ftjj^  g-iHot/figi(P)Rs/fig-i(R)Ps/ft_  (3.45) 
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Wavepackets  can  be  propagated  using  fewer  grid  points  in  the  nested  interaction 
picture  because  they  always  remain  centered  on  the  grid.  Translation  of  the  wavepacket  is 
corrected  for,  but  at  intermediate  points  in  the  calculations  the  wavepackets  still  may  be 
spread  by  the  effect  of  Uq.  Further  overhead  is  introduced  by  the  necessity  of  calculating 
the  expectation  values  (R)  and  (P)  for  every  time  step,  then  using  them  to  update  both  the 
state  vector  (3.43)  and  the  Hamiltonian  (3.45).  This  is  described  in  detail  in  Section  5.1.  In 
some  cases,  the  nested  Hamiltonian  may  be  simplified  using  what  Tannor’s  group  calls  the 
Heisenberg  api5roach.  The  Heisenberg  approach  requires  an  analytical  expression  for  the 
potential,  and  may  fall  victim  to  undersampling  of  the  potential.  The  underlying  reasoning 
and  some  computational  results  of  the  Heisenberg  approach  are  described  in  Section  4.2. 
The  most  accurate  and  reliable  version  of  the  nested  interaction  picture  is  the  simplest,  so- 
called  “sequential,”  approach.  The  sequential  method  accepts  the  computational  overhead 
described  above,  leading  to  computational  times  which  may  be  several  times  longer  than 
an  equivalent  Schrodinger-picture  propagation  in  one  dimension.  The  sequential  nested 
interaction  picture,  however,  does  enable  accurate  computations  using  reduced  grid  sizes 
compared  to  propagation  in  the  SchrOdinger  picture.  This  is  demonstrated  in  Chapters  V 
and  VI. 
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IV.  Alternative  Computational  Implementations  of  the  Interaction  Picture 

Three  approaches  to  the  interaction  picture  were  investigated  in  the  course  of  this  project 
that  failed,  either  to  yield  reliably  accurate  results,  or  to  promise  any  computational  advam 
tage  over  the  Schrodinger  picture.  The  sequential  approach  is  the  simplest,  and  is  the  equal 
of  the  Schrodinger  split-operator  method  in  stability  and  accuracy,  but  requires  the  same 
grid  size  and  more  computational  time  than  the  split-operator  method.  Tannor’s  so-called 
“Heisenberg”  approach  can  reduce  grid  and  computation-time  requirements,  but  requires 
an  analytic  expression  for  the  potential  and  is  not  accurate  for  many  interaction  poten¬ 
tials.  A  new  approach  was  also  developed  and  investigated,  and  named  the  “finite-basis 
approach.”  This  approach  proved  to  be  computationally  unstable. 

4.1  The  Sequential  Method 

The  non-nested  sequential  approach  was  investigated  early,  not  on  the  basis  of  any 
possibility  of  computational  advantage,  but  simply  to  verify  the  claim  of  Section  3.5.2  that 
Mpller  states  can  be  computed  by  single  interaction-picture  propagations.  The  sequential 
process  uses  the  Lanczos  propagation  scheme  of  Section  3,2  to  diagonalize  a  finite-basis 
approximation  of  the  interaction-picture  Hamiltonian, 

H/(f)  =  ut(f)HiUo(t).  (4.1) 

As  noted  in  Section  3.5.3,  the  diagonalization  process  for  this  Hamiltonian  requires  the 
same  grid  needed  for  the  Schrodinger  picture.  However,  it  is  comparatively  simple  to 
implement,  particularly  to  first  order,  and  serves  to  show  that  accurate  Mpller  states  do 
result  from  propagations  of  asymptotic  states  in  the  interaction  picture. 

This  is  demonstrated  by  the  computation  of  M0ller  states  of  the  square- well  potential 
of  Figure  4.1.  The  figure  shows  the  positive-momentum  reactant  Mpller  state  (x 
computed  in  the  Schrodinger  picture  using  a  split-operator  propagator,  based  on  the  pa¬ 
rameters  set  out  in  Table  4.1.  Figure  4.2  shows  the  same  Mpller  state  of  the  square  well, 
also  computed  using  the  parameters  stated  in  Table  4.1  in  the  interaction  pictme  with  a 
first-order  Lanczos  propagator  and  a  four-vector  Krylov  basis.  The  reactant  and  product 
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Figure  4.1  The  reactant  M0ller  state  4'+  (a;),  computed  in  the  Schrodinger  picture  for 
the  square-well  potential  represented  by  the  heavy  solid  line.  The  dotted  line 
represents  the  real  part  of  the  wavefimction;  the  thin  solid  line  its  modulus, 


X  {a.u.} 


Figure  4.2  The  reactant  M0ller  state  (x),  computed  in  the  interaction  picture  for 
the  square-well  potential  represented  by  the  heavy  solid  line.  The  dotted  line 
represents  the  real  part  of  the  wavefunction;  the  thin  solid  line  its  modulus, 

states  computed  in  the  two  pictures  can  be  compared  as  to  amplitude  and  phase,  according 
to  the  following  criteria.  After  Tannor  et  al.(49),  the  amplitude  error  is  given  by 


€a  —  1  —  1  5 


(4.2) 


and  the  phase  error  by 


(4.3) 
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Quantity 

Variable 

Value 

Coordinate  grid  spacing 

Ax 

.005 

Grid  size 

N 

4096 

Asymptotic  time 

r 

0.45 

Reduced  mass 

ij- 

1.0 

Time  step  size 

At 

0.001 

Initial  state  position 

Xo 

0.0 

Initial  state  momentum 

ko 

11.3 

Initial  Gaussian  width  parameter 

a 

0.25 

Table  4.1  Parameters  used  to  generate  the  M0ller  states  in  Figures  4.1  and  4.2.  All 
quantities  are  in  atomic  units. 


where  a  and  b  are  the  real  and  imaginary  parts  of  the  overlap  of  the  respective  wavefunctions- 
i.e.,  {'ipf  |‘0|)  =  a  -f  ib.  The  pairs  of  M0ller  states  illustrated  have  amplitude  error  ca  = 
3.7  •  10"^  and  phase  error  =  4.2  •  10~^. 

In  the  Schrodinger  picture,  the  process  of  computing  the  Mpller  states  involves  prop¬ 
agating  the  reactant  state  back  to  time  -r  and  the  product  state  forward  to  time  +r . 
As  illustrated  in  Figure  4.3,  this  drives  the  requirement  for  a  grid  larger  than  that  which 
would  be  needed  to  support  only  the  initial  and  M0ller  states. 

In  the  interaction  picture,  in  contrast,  the  intermediate  states  are  identical  to  the 
initial  states,  since  this  propagation  is  done  using  a  free-particle  potential.  The  M0ller 
states  evolve  directly  as  single  propagations  from  time  t  =  ±t  to  t  =  0.  Were  it  not  for 
the  internal  workings  of  the  sequential  method,  the  interaction  picture  would  thus  support 
the  entire  computation  on  a  smaller  grid.  However,  the  sequential  wavefunction  \'ip{t))  is 
computed  directly  as 


Hit))  I  =  vl{t)\i^{t))s 

=  vlit)v{t,o)Voio)Hmi 
=  vlit)V{t,o)H{o))s 

=  vUt)Hit))s-  (4-4) 
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{  n-e}  (x)A 


Figure  4.3  The  intermediate  state  {x 


g-i/fo(-T)|  ^  a.t  time  -r  = 


4.5  atomic  units. 


<x|exp{/H„t}|'P, 


Thus  at  time  ±t,  the  intermediate  states  derived  from  the  very  states 

iT))g  depicted  in  Figure  4.3,  and  therefore  require  the  same  grid  size  as  the  SchrSdinger 
picture. 

4-2  The  “Heisenberg”  Method 

Tannor’s  “Heisenberg”  approach  succeeds  in  bypassing  the  Schrodinger  picture  and 
trimming  the  computational  grid  requirements  for  Mpller  states  in  the  interaction  pic- 
ture(49).  Taimor  named  the  method  “Heisenberg”  because  of  the  similarity  of  the  equa¬ 
tions  of  motion  of  its  position  and  momentum  operators  to  the  Heisenberg  picture’s.  It 
is  not,  however,  a  true  Heisenberg  picture;  therefore  here  the  name  appears  in  quotation 
marks.  The  “Heisenberg”  approach  is  based  on  manipulation  of  an  analytic  expression 
of  the  Hamiltonian,  using  a  number  of  operator  identities.  The  construction  of  the  time- 
evolution  operator  in  the  method  begins  with  a  reformulated  potential  and  Hamiltonian. 

4-2.1  The  “Heisenberg”  Potential.  Consider  a  full,  iV-dimensional,  charmel 
Hamiltonian  in  Jacobi  coordinates, 

j  =  l  W 

where  the  operator  R5  includes  both  internal  and  external  coordinates.  Let  a  generic 
operator  Hq  be  defined  as  the  sum  of  any  portion  of  the  kinetic-energy  operator  and  some 
analytic  portion,  $  (Rs)  =  J2n=i  ^jri^sj’  potential  operator: 

H3=E^+®(Rs).  (4.6) 

The  equation  of  motion  for  an  interaction-picture  operator  Oj  (t)  =  e»Hot//iQ^g-iHo«/A  jg 

ihj^Oi  =  ih(^^iio'^Oi  +  ihOi(^^Jio 

=  [0/,H2].  (4.7) 
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Therefore,  the  time  evolution  of  a  position  operator  R/,„  =  ***“*/^ 

by 


^  [H2,R/m] 


h 


L  ■p2 


h  2/^i  . 


IpiHot/ft 


L  p2 
^Sj 


J^^+^(ns),Rsr, 

1=1  J 


2h 


i=i 


p2 

-^,Rsm 


J_  iHot/n 

2h 

_L  iHot/fi 

2h 

1  „ 

P  /m* 


[Pk, 

L  Mm 

/z^ 

\  Mm 


RSn 


P5m  i 


The  derivation  of  equation  (4.8)  uses  the  commutator  relations 


[P j  ?  Pm]  —  (P j  P j  Pm  PmPj  Pj  )  ^jm 

=  (PmPmPm  “  PmPmPm  "1“  PmPmPm  ””  PmPmPm) 
=  (P m  [Pm?  Pm]  [Pm 5  Pm]  Pm) 

=  —i2flPrn^ 


[$  (Pj)  ^  Pm]  —  Oj 


and 


is  given 


(4.8) 

(4.9) 

(4.10) 

(4.11) 
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The  time,  evolution  of  a  momentum  operator,  'Pim  =  is 


dt^ 


h 


y  pi  +  $  (R5)  , 

AT  M 

j=l  n=l 
M 

n=l 

M 

=  [R;j,5,P5m] 

n=l 

M 

n—1 

=  (RmS)  6-*““*/'' . 


(4.12) 


If  the  constraint  V$m  (Rms)  =  0  is  imposed,  then  we  have 

s'*"”"®- 


(4.13) 


Equation  (4.13)  indicates  that  each  component  P/m  of  the  momentum  operator  P/  is 
constant  in  time,  so  that  P/  {t)  =  P/  (0)  =  P5.  Equation  (4.8)  can  thus  be  integrated 
directly,  to  yield 


R'/m  (^) 


(4.14) 


Hence,  by  induction, 


R?™(<)  = 


Rsm  H - PSti 

/^m 


(4.15) 
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for  any  natural  number  n,  and 


V  (R/^)  =  v(Rsm  +  —Psm)  (4.16) 

for  any  analytic  potential  function  ^(R). 

4.2.2  The  “Heisenberg”  Hamiltonian  and  Evolution  Operator.  Using  the  reason¬ 
ing  of  section  4.2.1,  if  the  coordinate  operator  is 

Ri{t)  =  Rs  +  —Ps,  (4.17) 

then,  after  two  applications  of  the  operator  identity 


gARe-A  =  B  +  [A,  B]  +  I  [A,  [A,  B]]  +  ^  [A,  [A,  [A,  B]]]  +  . . . ,  (4.18) 


the  coordinate  operator  in  the  PR-adapted  nested  interaction  picture  (Section  3.5.4),  is 


^  g-i(Rs)Ps/ftg-t(Ps>Rs/^  ^i{Fs)ns/ll^i{Rs)Ps/H 

^  g-i<Rs>Ps/^  /^Rs  +  —  + 

V  Ai  A^  / 

=  R5  +  ^  +  i^  +  (R). 

/X  /i 

Therefore,  the  PR-adapted  Hamiltonian  for  a  one-dimensional  analytic  potential  is 


(4.19) 


H"  =  iijei(P)P-s/h^i{R)Ps/fi 

=  hJrs  +  —  +  {R)  +  ^).  (4.20) 

The  time-evolution  operator  is  now  approximated  using  the  iterated  Lanczos  method. 
Krylov  basis  vectors  are  built  according  to  equation  (3.6),  by  repeated  operation  of  the 
Hamiltonian  H"  on  the  wavepacket  t/'  (t).  The  computational  savings  of  the  “Heisenberg” 
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approach  compared  to  a  sequential  approach  vary  depending  on  the  chosen  partition  of  H 
into  Ho  and  Hi. 

j^,2.3  A  Reformulated  ‘^Heisenberg’’  Hamiltonian.  The  manipulations  of  Section 
4.2.1  can  lead  to  simpler  formulation  of  the  evolution  and  Mpller  operators  in  certain  cases. 
In  constructing  a  particular  interaction  picture,  the  choice  of  the  operator  Hq  is  free. 
However,  that  choice  immediately  determines  the  remaining  portion  of  the  Hamiltonian  to 
be 


H^  -  H^  -  H^ 

=  E  ^  + 

j=L+l 

=  E  +  (4.21) 

The  full  interaction-picture  Hamiltonian  is  then  written, 


H] 


=  e 


N  p2 

iHot/h  (  _i_  T/' 


E  if 


U=L+1 


.  ^  P2  . 

^  giHot/ft  j  ^  :L£i  j  g-iHot/fi  ^Hotlhyi  g-iHot/n 

vj=L+l  "^^3 


N 


=  V  V'  (R/) 

j=L+l  ^^3 


N.  ^  /  i  \ 

=  y  -mot/h^Hot/hp  -mot/n  +  V  f  Rsm  +  —^Sm  ] 

JV  M  ,  . 

j=L+l  ^^3  m=\  \ 

JV  M  f  . 

=  E  +  + 

j=I,+l  ^^^3  m=:l 
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The  first  term  of  the  expansion  of  the  time-evolution  operator  is  therefore 


U7(t,to)  =  exp 


ij  («')*'} 


(4.23) 


'  r  ^  1  M  /  , 

=  exp  I  -!j—  ^  C  ('R'Sm  (to)  + 

"  j=L+l  m=l  ^ 


This  form  of  the  operator  is  amenable  to  split-operator  propagation  for  those  potentials 
V'  for  which  an  analytic  expression  is  available.  Exponential  potentials  such  as  Morse 
oscillators  are  particularly  tractable  using  this  approach.  Let  us  now  examine  the  results 
of  certain  choices  of  interaction  picture. 


4-2.4  Examples  of  “Heisenberg”  Interaction  Pictures. 


4.2.4. 1  The  One-Dimensional  Case.  In  one  dimension,  the  Hamiltonian  is 


given  by 


h='L^  +  v(o. 


(4.24) 


Ho  =  H„  = 


2m  ’ 


(4.25) 


Hi=  V(r). 


(4.26) 


The  Moller  operators  are  time-evolution  operators,  as  shown  in  Section  3.5.2.  For 
small  time  intervals,  the  time-evolution  operator  is  approximated  by 


Vj{t,to)  =  exp|-^jf’*H/(t')dt'| 
exp  /  (ts+ 


ks  . 


(4.27) 
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This  form  is  fast  and  simple  to  apply  directly,  but  again  requires  an  analytic  expression 
for  the  potential  V  (x). 


Two-Dimensional  Cases.  Suppose  we  have  a  two-dimensional 
Hamiltonian  in  Jacobi  coordinates, 


^  +  y^lr^,  a,). 


(4.28) 


As  mentioned  earlier,  a  number  of  different  partitions  of  the  Hamiltonian  into  the  portions 
H2  and  are  available  in  systems  with  two  or  more  degrees  of  freedom. 


•  Consider  the  choice, 


(4.29) 


H]  =  Ily).  (4.30) 

In  this  interaction  picture,  the  Mpller  operators  are  phase-shifted  time  evolutions. 
The  short-time  evolution  operator  is 

=  expl^jf  Hj(t')df'| 

«  exp|^[^k|i^-FK,(r.,  +  ek,,,R,  +  0k^)  (4.31) 

where 


h  [tp  +  (t  ~  tp)  /2] 


(4.32) 


This  evolution  operator  can  be  applied  in  split-operator  fashion,  provided  always 
that  the  interaction  potential  is  available  in  analytic  form. 
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•  Suppose  instead  we  choose, 


2m^  2/i^ 


(4.33) 


H'l=Y^{r^,Ry).  (4.34) 

In  this  interaction  picture,  the  short-time  evolution  operator  is  approximately 

=  expj-^^  Hj(i')d^'} 

exp  (*‘7'S  Ro'S  +  ^k/j^s)  I ,  (4.35) 

(where  6  is  given  once  again  by  equation  4.32)  which  can  be  applied  directly  in  the 
coordinate  representation,  given  an  analytic  form  of  the  potential. 
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The  M0ller  state,  given  by 


1^+)/  =  k+>5 

=  lim  V\{t)Va{t)\^PJs 

<— )■— oc  ' 

=  lim  (t)  (t) 

t— »-oo 

=  lim  g-iH-yt/figiHjt/ft  1^  J 
=  lim 

t—yoo 

«  lim  e-»H2‘/2'ie-»2H?t/2fig-iH2t/2figiH2:t/2figi2V7„,t/2ftgiH2t/2ft 

t— >oo 

=  lim  e-»2H? (|)i/2ftgi2V7„,^(|)t/2ft 
t— >oo 

=  lime-2H?Wt/fte*2V7n.,W‘/ft|^.j^ 
t—>oo 

r/At  r/At 

«  n  n  iV'in)/ 

n=l  m=l 

-i2At  ^  \  'i 


=  n^’^pj— 

n=l  ^ 


V,  Ts  +  '-^ks 


n  ®^p  { 1 ’ 


(4.36) 


(employing  the  Zassenhaus  identity  (3.2)  involves  two  propagations  in  this  case,  since 
Ho  #  Ha.  The  final  line  of  equation  (4.36)  is  only  of  use  for  analytic  forms  of  the  potentials 
Vy  and  It  should  also  be  noted  that  the  allowable  size  of  the  time  step  At  is  halved 
in  this  case  because  the  Hamiltonian  is  doubled  in  the  propagation  operator. 


•  A  final  interesting  partition,  choosing  Hq  as  the  full  Hamiltonian,  results  in  a  true 
Heisenberg  picture.  In  this  special  case  of  the  interaction  picture,  state  vectors  are 
time  independent; 

\^)h  =  1^  (t))5 

=  giHt/fig-iHt/fi  1^  (0))^ 

=  |V’(0))5-  (4-37) 
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Therefore,  time-dependent  propagation  methods  do  not  apply  in  the  way  discussed 
thus  far.  Instead,  operators  must  be  evaluated  directly  as  functions  of  time. 


Oh  (t)  = 


(4.38) 


The  Mpller  operators  are  only  defined  at  t  =  0,  so  they  are  the  same  as  in  the 
Schrodinger  picture. 


(4.39) 


The  problem  of  computing  Mpller  states  is  thus  formally  identical  in  the  Schrodinger 
and  Heisenberg  pictures.  The  Moller  operator  could  be  decomposed  into  its  two 
evolution-operator  components,  and  each  operation  performed  in  the  Heisenberg  pic- 
tmre  using  a  variant  of  the  Lanczos  technique,  but  this  approach  has  no  apparent 
advantage  over  a  judiciously  chosen  interaction-picture  propagation. 


4-£-5  Computational  Problems  of  the  “Heisenberg”  Approach.  First-  and  second- 
order  Lanczos  propagators  were  developed  to  demonstrate  the  “Heisenberg”  interaction 
picture.  The  results  were  good  for  the  exponential  potential  functions  V  (x)  = 
already  demonstrated  by  Tannor’s  group(49),  but  not  at  all  reliable  for  realistic  interaction 
potentials.  The  problem  may  lie  in  the  fact  that  the  effect  of  the  potential  is  evaluated 
entirely  in  the  Krylov  subspace,  using  only  a  few  sampling  points  instead  of  the  many  used 
when  the  potential  is  sampled  at  each  grid  point  in  coordinate  or  momentum  space.  In 
a  one-dimensional  PR-adapted  nested  “Heisenberg”  interaction  picture,  for  example,  the 
Hamiltonian  is  evaluated  as 


where  the  argument  y  (i)  =  x  +  ^  +  (x)  +  is  diagonalized,  leading  to  the  potential 
being  evaluated  at  the  Ritz  values  yi  (t)  =  {qi  |y  (t)  |gi)  Figure  4.4  shows  the  dispersion  of 


Figure  4.4  The  coordinate  representation  of  the  first  seventeen  Ritz  values  of  a  Gaus¬ 
sian  wavepacket  propagated  in  the  “Heisenberg”  interaction  picture  using  a 
twenty-vector  Krylov  basis,  with  a  Gaussian  well  potential.  The  progression 
of  the  Ritz  values  over  time  is  represented  with  open  circles.  The  Hamilto¬ 
nian  is  evaluated  as  the  value  of  the  potential  V  (y),  represented  by  the  solid 
curve,  where  only  the  Ritz  values  of  y  are  used  instead  of  the  relatively  much 
more  densely  distributed  coordinate  values  x. 

the  Ritz  values  for  a  localized  scattering  potential  as  propagation  time  progresses.  Since 
the  Hamiltonian  is  approximated  by  evaluating  the  potential  fimction  at  the  Ritz  values, 
the  Ritz  values  are  plotted  in  conjunction  with  V  (x).  As  time  progresses,  some  Ritz  values 
follow  the  dispersing  wavepacket,  while  some  remain  in  the  region  where  the  potential  is 
significant.  This  opens  the  computation  to  loss  of  accuracy  as  the  resulting  basis  may 
come  to  bear  less  and  less  of  a  relationship  to  either  the  Hamiltonian  or  the  wavefunction. 
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This  is  a  speculative  explanation  for  the  observed  unreliability  of  the  technique.  Of  several 
analytic  potentials  examined  using  this  technique,  only  an  exponential  potential  yielded 
accurate  propagated  wavepackets  using  the  “Heisenberg”  technique. 


4-3  The  Finite-Basis  Approach 

The  third  approach  tried  and  abandoned  was  named  the  finite-basis  approach.  It  had 
the  advantage  of  speed,  based  on  the  virtual  elimination  of  the  requirement  to  perform 
FFTs  by  performing  the  entire  computation  in  a  Krylov  subspace  based  on  either  the 
coordinate  or  momentum  representation  of  the  wavefunction,  rather  than  both.  However, 
this  approach  could  not  be  made  computationally  stable. 


4.3.1  Derivation.  In  the  iterative  Lanczos  method(44),  the  evolution  operator 
U/  (t)  is  expressed  in  a  finite  basis  of  dimension  M  N,  where  N  is  the  grid  size(49-51). 
In  matrix  notation. 


{xi  lip  {t  +  At)) 


{xn  \tp  {t  -h  At)) 


{«!  \qi)  (a;i  Iqm) 

•  ■■■  ;  (4.41) 

{xN  kl)  •  •  •  {XN  Iqm) 

'  (gi|U(t  +  At,t)ki}  (9i|U(t  +  A<,t)kM) 

{qM\V  {t  +  At,t)\qi)  ■■■  {qM\lJ{t  +  ^t,t)\qM) 


{qi  ki) 

X  ; 

{qM  laJi) 


{qi  kiv) 

{qM 


{xi  \ip  (t)) 


{XN  IV’  (<)) 
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The  M  basis  vectors,  |9„),  in  the  Lanczos  method,  are  derived  from  orthonormalized 
Krylov  vectors 

\<j>^)  =  H?  \xl,  (t))  =  Uj  (t)  V'^Uo  (t)  1^  {t)) .  (4.42) 

This  expression,  if  computed  directly  as  in  the  sequential  method,  involves  two  free-particle 
type  propagations  over  the  full  time  t,  not  just  the  incremental  time  At:  first  by  Uq  (t)  and 
then  by  Uq  (t).  This  double  propagation  must  be  done  M  times  for  each  time  step  in  order 
to  obtain  all  of  the  Krylov  vectors.  Furthermore,  the  same  grid  size  is  required  as  in  the 
Schrodinger  picture,  since  the  full  spreading  effect  of  Uq  (t)  is  ultimately  incurred  before 
being  reversed  by  (t).  This  direct  approach  for  computing  the  Krylov  vectors  must 
therefore  be  avoided,  if  propagation  in  the  interaction  picture  is  to  have  any  advantage 
over  the  Schrodinger  picture. 

The  need  to  propagate  for  the  full  time  t  and  back  again  can  be  finessed  by  performing 
the  propagations  using  Uq  (±At)  instead  of  Uq  (it).  The  interaction  Hamiltonian  at  any 
time  tj  ==  tj-i  +  At  is  related  to  the  previous  time’s  Hamiltonian  as 

Kj{tj)  =  Vl  (tj)  Wo  (tj) 

_  giHo  (tj  ~  (tj  _  i  -\-At)/h 

_  ^iHoAt/h  ^—iHoAt/h 

=  (tj-i)  e-***'>^*/*.  (4.43) 

The  previous  Hamiltonian,  H/  (tj-i),  has  already  been  diagonalized  in  its  own  finite  basis, 
and  may  be  transformed  to  the  full  momentum  representation  by  means  of  a  non-square 
matrix  Q,  whose  elements  Qnm  =  (kn  are  computed  as  part  of  the  Lanczos  tri- 

diagonalization  process.  This  enables  H/  (tj)  to  be  represented  in  turn  in  an  M-dimensional 
subspace,  and  then  applied  to  a  wavepacket  in  the  momentum  representation. 

For  each  time  step  tj,  the  first  Krylov  vector  is  chosen  to  be 

9j)  =  IV>(tj-l)),  (4.44) 
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and  is  computed  in  its  momentum  representation,  (k  To  obtain  each  subsequent 


basis  vector, 


it  is  necessary  to  calculate  the  vector 


#i_i)  =  H,((i)|4_ 


(4.45) 


In  the  momentum  representation. 


U,) 


MM 

=  /  dk'  (k  k')  £  (k'  kr*)  £  (ej'  |H,  ((^_,)| 

J  ^  ^  n'=l  n"=l 

X  j  dk"{^e^^\k")  J  dk'"  (fc"  k'")  {k'"  .  (4.46) 

The  linearly  independent  Krylov  basis  vectors  |^)  are  orthonormalized  to  obtain  the 
basis  vectors  Thus,  the  basis  vectors  are  calculated  using  propagations  over  just  the 
durations  ±At,  enabling  a  reduction  in  grid  size  and  eliminating  FFTs  entirely  from  the 
propagation  process. 

Further  examination  of  the  process  for  calculating  subsequent  basis  vectors  reveals 
another  benefit  of  pairing  the  interaction  picture  with  the  Lanczos  propagation  method. 
The  vector  is  generated  according  to  equation  (3.9), 


|9j> 

II 

(H/(t 

_  1  J 

r  1  . 

1 

^n— 

=  {tj)  [H/  (t,)  -  an-2  \(i-2)  -  0n-2 

«n-l  \4i-l)  -  ^n-1  1^-2)! 


(4.47) 


A 


^  [Hj  {tj)  Hj  (t,)  \qi_,)  -  a„_2H/  (t,)  \qi_^)  -  (tj)  |q^„_3)] 

nPn—1 


1 

■77 


/?n 
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The  vectors  ),  for  1  <  m  <  M,  are  given  by 


I  A)  =  »,(«,)!<) 

=  H/  (tj)  |H/  (tj)  —  Oim-l  Qm-l)  ~  ^m-1  } 

=  ^{Hi{tj)\ji,)-arn-ili,.,)-l3m-i\pL-2)}-  (4-48) 

Pm  ^ 

Then  for  1  <  n  <  M,  equation  (4.47)  becomes 


=  I  -  Q^n-l  kn-l  (ij))/ “ /^ra-l  |9?i-2  »  (4.49) 

where  the  values  of  i^n-2)  known  from  previous  calculations.  Only  one 

new  vector,  H/  (tj)  |p^_2^)  is  computed  for  the  first  time.  The  sum 

\li-l)  =  {H/  (*i)  \pL-2)  -  «n-2  \li-2)  -  ^n-2  ll^-s)}  (4-50) 

is  saved  for  future  use. 

It  is  now  clear  that  the  operation  of  the  Hamiltonian  is  going  to  be  evaluated  only  in 
expressions  of  the  form  of  the  first  term  of  equation  (4.50).  Such  vectors  can  be  evaluated 
as  outlined  above,  by  recoiurse  to  the  previously  diagonalized  operator  H/  (tj-i): 

H/(tj)^_2)  =  ^{H/(t,)H/(t,)|?4_3)-a„-3H/(t,)^_3) 

-  ^n-aH/  (tj)  |j4-4)} 

=  -j—  {tj-i)  e-iHoAt/ftgiHoAt/fijj^  g-iHoAt/fi  pJ  _  \ 

0n-2  *' 

-  (t,_i)  ?4_4)}  (4.51) 

=  f  U\  _  . 

Pn  ^ 
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The  terms  in  equation  (4.51)  represent  the  vectors 

(4.52) 

and  for  1  <  m  <  M, 

14,)  =  H/  {tj-i)  14,)  ,  (4.53) 

which  turn  out  to  be  a  more  convenient  form  in  which  to  store  the  information  contained  in 
the  vectors  taking  advantage  of  the  cancellation  of  evolution  operators  which  occurs 
in  repeated  operations  of  H/(tj_i).  The  operations  contained  in  the  braces  in  equation 
(4.51)  are  all  performed  in  the  Krylov  subspace  of  the  time  t  =  tj-i,  and  the  result 
transformed  to  the  momentmn  representation  and  propagated  by  —At.  The  propagated 
vector  is  then  added  to  the  remaining  terms  of  equations  (4.49)  and  (4.50)  in  momentum 
space.  Only  M  + 1  FFTs  are  required  over  the  entire  propagation  to  compute  the  first  M 
basis  vectors  if  Ho  and  Hi  are  respectively  functions  of  position  and  momentum  alone,  as 
opposed  to  2M  FFTs  per  time  step  required  by  the  Lanczos  method  in  the  Schrodinger 
picture. 

4-3.2  Implementation  Problems  and  Possible  Solutions.  The  finite-basis  algo¬ 
rithm  was  implemented  in  a  one-dimensional  propagator,  using  a  first-order  truncation  of 
the  time-evolution  operator  and  no  re-orthogonalization  of  the  basis  vectors.  In  this  crude 
form,  the  method  falls  victim  to  two  idiosyncracies  of  the  iterative  Lanczos  approach. 

First,  equation  4.46  implicitly  assumes  that  the  4^  span  the  full  IV-dimensional 
space.  This  approximation  that  \i)  (tj+i))  can  be  expressed  as  a  linear  combination  of  the 
1 4^  is  adequate  for  a  single  time  step,  but  its  error,  compounded  over  several  time  steps, 
is  sufficient  to  cause  the  pure  finite-basis  method  described  above  to  become  unstable  for 
long-time  propagations  of  the  type  required  for  reactive  scattering  calculations.  The  loss  of 
accuracy  begins  immediately  with  the  highest-order  Krylov  basis  vectors,  so  increasing  the 
basis  size  delays  the  onset  of  inaccuracy  for  the  wavefimction  itself,  but  does  not  prevent 
it  permanently.  Once  the  degradation  progresses  far  enough,  too  few  appropriate  basis 
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vectors  remain  to  represent  the  wavefunction  accurately.  Figure  4.5  shows  the  onset  of 
degradation  in  wavefunctions  constructed  in  Krylov  bases  of  various  order. 


Time  {  atomic  units  ) 


Figure  4.5  Increasing  error  with  propagation  time  for  Gaussian  wavepackets  on  a  linear 
potential,  using  a  first-order  finite-basis  nested  interaction  picture  propagator 
with  various  Krylov  basis  dimensions  M.  The  reference  wavefunction  is 
computed  using  the  sequential  non-nested  interaction  picture  described  in 
Section  4.1. 

Second,  the  Lanczos  tri-diagonalization  algorithm  suffers  from  a  tendency  for  the 
orthogonality  of  the  basis  vectors  to  degrade  with  repeated  iterations.  This  results  from 
a  well-known  numerical  problem  of  Gram-Schmidt  orthogonalization,  which  can  be  dealt 
with  when  necessary  by  periodically  re-orthogonalizing  all  or  some  of  the  basis  vectors(72— 
75). 


Three  possibilities  for  stabilizing  the  finite-basis  algorithm  were  considered: 

•  Periodic  “resetting”  of  the  approximation  by  performing  a  direct-propagation  step 
using  FFTs  on  the  momentum  grid.  This  lengthens  the  lifetime  of  the  propagation, 
but  ultimately  the  propagation  still  fails  catastrophically  (Figure  4.6).  Degradation 
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Figure  4.6  Degradation  in  accuracy  of  propagated  wavepacket  with  time.  Propagations 
were  done  on  a  linear  ramp  potential,  using  the  finite-basis  algorithm,  “re¬ 
set”  periodically  by  performing  a  time  step  using  the  sequential  non-nested 
interaction  picture.  The  reference  wavefunction  ^ref  is  computed  entirely  in 
the  sequential  non-nested  interaction  picture. 
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is  not  entirely  eliminated  even  when  every  other  time  step  is  performed  using  the 
inefficient  non-nested  sequential  method,  as  reflected  in  the  gradual  destruction  of 
the  state  vector  illustrated  in  Figure  4.7. 


Figmre  4.7  Even  alternating  every  other  iteration  between  nested  finite-basis  and  non¬ 
nested  sequential  interaction  pictures  results  in  a  gradual  but  continuous 
degradation  in  accuracy  of  the  propagated  wavepacket  as  compared  to  '^ref, 
obtained  entirely  by  using  the  sequential  non-nested  propagator. 

•  Re-orthogonalization  of  the  basis  vectors  to  allow  larger,  more  stable  basis  sets. 
This  was  attempted,  and  failed  to  have  any  effect  on  the  degradation  of  accuracy  in 
propagated  wavepackets. 

•  Implementation  of  higher-order  Magnus  representations  of  the  time  evolution  opera¬ 
tor  to  improve  the  accuracy  of  transitions  from  one  time  step’s  basis  to  the  next,  and 
allow  longer  propagations  with  less  error.  This  would  only  postpone  the  degradation 
by  allowing  longer  propagations  with  the  same  number  of  time  steps,  so  second-order 
Lanczos  propagation  was  never  applied  to  the  finite-basis  approach. 
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V.  Software  Implementation  and  Validation 
5.1  Code  Implementation 

The  core  activity  of  this  project,  described  in  Chapter  VI,  involves  computing  S- 
matrix  elements  by  generating  Mpller  states  in  the  interaction  picture,  then  evaluating 
the  Fourier  transform  of  the  time-dependent  correlation  function  (2.24),  computed  in  the 
Schrbdinger  picture  with  absorbing  boundary  conditions.  The  Mpller  states  are  generated 
using  a  Lanczos  propagator  based  on  the  second-order  Magnus  expansion  of  the  time- 
evolution  operator  to  propagate  product  and  reactant  wavepackets  from  t  —  ±t  to  t  =  0. 
The  Tannor  sequential  PR-adapted  nested  interaction  picture,  introduced  in  Section  3.5.4, 
is  used.  For  comparison,  Mpller  states  are  also  computed  in  the  Schrodinger  picture,  using 
the  split-operator  method  of  Section  3.1.  The  SchrSdinger  spilt-operator  method  is  faster 
and  more  commonly  used  than  the  Lanczos  method  in  the  Schrodinger  picture,  and  thus 
presents  the  primary  benchmark  for  comparison  of  computational  effort  required  to  obtain 
similar  results  in  the  interaction  picture. 

All  code  for  calculating  the  Mpller  states  is  written  in  C-|— 1-.  Several  legacy  subrou¬ 
tines  in  standard  C  are  employed,  including  Basic  Linear  Algebra  Subprograms  (BLAS) 
and  CLAPACK  linear-algebra  subroutines  from  Netlib,  and  a  matrix  diagonalizer  and  FFT 
from  Numerical  Recipes  in  C  (76-79). 

The  Mpller  states  are  passed  to  a  modified  version  of  a  Fortran  program,  developed 
by  R.  S.  Calfas  to  generate  correlation  functions,  using  split-operator  propagation  with 
absorbing  boundary  conditions  in  the  Schrodinger  picture(23,80).  S-matrix  elements  are 
then  derived  from  the  correlation  functions  using  another  C-t— I-  application.  The  Delta 
C-b-l-  compiler  and  Silicon  Graphics  Fortran  compilers  were  used  on  MIPS  3000  and  higher- 
performance  Silicon  Graphics  workstations. 

5.1.1  Initial  Conditions.  The  interaction-picture  Mpller  states  are  computed  in 
the  PR-adapted  interaction  picture(49)  defined  by 


The  channel-packet  method  begins  with  four  initial  states, 


^  in  !  out 
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(5.2) 


made  up  of  one  reactant  (in)  and  one  product  (out)  state  in  each  of  the  two  channels  (to  the 
left  (A)  and  right  (p)  of  the  interaction  region  of  the  potential.  These  initial  wavepackets 
are  constructed  analytically  in  the  Schrodinger  picture  as  Gaussians 

(5.3) 

(the  coordinate  representation  of  Equation  2.5)  centered  near  the  interaction  region  at 
i  =  0,  but  are  the  same  as  the  corresponding  interaction-picture  states  (x  (±'^))^ 

at  time  -hr  for  products  and  — r  for  reactants.  The  corresponding  PR-adapted  intermediate 
states, 
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(5.4) 


are  then  calculated,  using  the  initial  conditions  {x)  —  xq  and  (k)  =  fco*  Subsequently, 
new  values  of  (x)  and  {k)  must  be  computed  after  each  time  step  and  used  to  update  the 
evolving  wavepacket.  At  t  =  to?  the  end  of  the  channel-packet  method’s  first  propagation, 
the  M0ller  states  are  converted  from  the  nested  interaction  picture  to  the  Schrodinger 
picture,  as 


(a:|^±)5 


(x 


g-i(k)Xgi{x)k 


(5.5) 


5,L2  Computation  of  M0ller  States,  The  short-time  evolution  operator  is  ap^ 
proximated  by  a  second-order  Magnus  expansion, 
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where  midpoint  quadrature  is  employed  to  approximate  the  first-order  term,  and  trape¬ 
zoidal  integration  for  the  second-order  term.  The  geometry  of  the  derivation  of  the  second- 
order  term, 

1  /‘tfc+i  ft' 

B  (4+1,4)  =  ^l  dt'j^  df  [H  (<')  ,H  (t")] ,  (5.7) 

is  illustrated  in  Figure  5.1.  Two-dimensional  trapezoidal  quadrature  operates  in  a  super- 


Figure  5.1  Calculation  of  the  two-dimensional  trapezoidal  integral. 

space  which  adds  the  time  dimensions  t'  and  t"  to  the  Hilbert  space  containing  the  com¬ 
mutator.  The  hypersurface  of  integration,  z  (t',  t")  =  [H  {t') ,  H  (t")] ,  is  approximated  by  a 
plane.  The  integral  is  the  volume  between  the  plane  and  the  triangle  containing  the 

points  {t'  =  t"  =  4, z  =  0),  (t'  =  t"  =  4+1, ^  =  0),  and  {t'  =  4+i, —  tk,z  =  [H  (4+i) , 
H  (4)]  =  C).  The  integral  is  easily  performed  in  the  coordinates  x  =  t"  —  tk,  y  =  4+i  ~ 
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where  6  =  At  =  tjc  —  tk-i  •  Thus, 


Jl_c 

(M! 
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[H(4+i),H(ife)]' 


(5.8) 


The  result  of  applying  the  evolution  operator  to  a  PR- adapted  wavepacket  (^x  j 


at  time  t  —  tj  is 


(a:)  =  (x  , 


(5.9) 


a  representation  of  the  wavepacket  at  time  t  =  tj+i,  but  in  the  interaction  picture  defined 
by  the  displacement-boost  operator  corresponding  to  time  t  —  tj.  The 

wayepacket  is  brought  into  the  correct  interaction  picture  by  first  updating  the  expectation 
values, 


(x),„  =  {x),  +  (P+‘|>'l«^'"'>.  (6'W) 


and 


(k),+i  =  (k),+<e+‘|k|c*’>,  (5.11) 
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then  computing, 


{x\rp^+^y'  = 

=  ei{^)j+ikg-mj+ix^mjx^-i{^)^k^j+i 

_  e*(£|x|C>  fegi<x>jfcg-i(€|k|?)  ®g-i(x)j.fc^j+l 

=  ei<C|x|e>  fcg-i(e|k|e)  Xgi(x),fcg-i(?|k|0  (x>,g-j(x>jfc^i+l  (3.) 

=  gi«|x|Ofee-®<^l‘‘lf>^e"*<^l‘‘l^^  (x).  (5.12) 


The  derivation  of  equation  (5.12)  uses  the  lemma, 

gakg6x  ^  gaxgfckgiat^ 


(5.13) 


which  is  a  corollary  of  the  Zassenhaus  formula  (3.2),  since 

gOx+fek  _  gaXg6kg-i[ax,6k]  _  g6k+ox  _  ^bk^aXgA  [ax,6k]  ^  (5.14) 


5.2  Validation  with  Square  Potentials 

To  confirm  the  capability  of  the  interaction-picture  approach  to  compute  accurate 
Mpller  states  and  S-matrix  elements,  the  process  can  be  tested  using  an  asymmetric  square- 
well  potential,  with  one  asymptotic  potential  energy  higher  than  the  other  (Figure  5.2). 
This  potential  is  chosen  because  its  transmission  and  reflection  coefficients  can  be  calcu¬ 
lated  analytically,  and  it  exercises  the  capability  of  the  channel-packet  method  to  deal  with 
disparate  asymptotic  energy  levels.  Square  potentials,  however,  cannot  be  represented  with 
complete  accuracy  on  a  discrete  grid.  No  finite  grid  can  support  the  true  vertical  slopes 
that  characterize  such  a  potential;  however,  a  successful  method  will  be  able  to  demonstrate 
convergence  toward  the  analytic  solution  as  the  potential  is  more  accurately  approximated. 

5.2.1  Sources  of  Numerical  Error.  All  time-dependent  propagation  techniques 
are  prone  to  error  resulting  from  the  discretization  of  continuous  events.  Position  and 
momentum  are  considered  as  finite  sets  of  non-contiguous  points,  usually  separated  by 
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Figure  5.2  The  asymmetric  square-well  potential  used  to  validate  the  propagators.  The 
well’s  discontinuities  are  at  ±1  atomic  unit,  the  well  bottom  is  at  -100  atomic 
units  of  energy,  and  the  left  and  right  asymptotes  are  at  0  and  50  atomic  units. 


constant  spacings  Aa:  and  Ak  to  accommodate  the  FFT.  Time  also  is  divided  into  discrete 
intervals  At.  Decreasing  the  size  of  any  of  these  intervals  increases  the  accuracy  of  the 
numerical  approximation,  while  increasing  the  numerical  effort  required.  Discontinuous 
potentials  such  as  the  square  well  pose  special  numerical  problems  of  their  own.  When 
expressed  on  a  discrete  coordinate  grid,  what  should  be  square  becomes  trapezoidal,  as  the 
infinitesimal  distance  across  the  discontinuity  is  stretched  out  to  a  finite  Ax.  This  effect 
is  illustrated  in  Figure  5.3. 


X  (a.u.) 


Figure  5.3  The  asymmetric  square  well  potential,  to  successively  better  approximations 
resulting  from  finer  grid  spacings. 

5.2.2  Analytic  Calculation  of  Transmission  and  Reflection  Coefficients.  The 
analytic  calculation  of  the  probabilities  of  transmission  and  reflection  is  a  straightforward 
but  interesting  exercise  based  on  continuity  of  wavefunctions  and  their  derivatives  across 
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the  discontinuities.  Bohm’s  lucid  treatment  of  the  simple  square  well  is  easily  adapted  to 
accommodate  one  asymptote  of  the  potential  at  nonzero  energy(81).  For  a  square  well 
with  left  asymptote  at  energy  =  0,  right  asymptote  at  E  =  Vi,  and  well  of  width  2a  at 
potential  —Vq,  the  transmission  coefficient  from  left  to  right  is 


Ti{E)  = 


+  +  2AB  cos  (4fc2a)) 


h 

fci’ 


(5.15) 


where 


A  = 


(5.16) 


(5.17) 


ki  = 


y/2^E 

h 


(5.18) 


A:i  = 


y/2^  {E  +  Vq) 


(5.19) 


and 


,  ^2fi{E-V,) 

h  • 


(5.20) 


The  probabihty  of  reflection  for  a  wavepacket  traveling  from  left  to  right  is  of  course, 


i?i  =  l-ri.  (5.21) 

The  analytic  transmission  and  reflection  coefficients  for  the  square  well  of  Figure  5.2  are 
shown  in  Figures  5.4  and  5.5. 
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Figure  5.5  Reflection  coefficient  for  the  asymmetric  square  well. 
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5.2.3  Numerical  Approximations  of  S-Matrix  Elements.  The  square  and  trape¬ 
zoidal  potentials  dealt  with  in  this  chapter  are  formed  along  the  lines  of  Figure  5.3,  with 
anchors  at  the  leftmost  points  of  the  well  and  right  asymptote.  Square  potentials  con¬ 
structed  in  this  manner  approximate  the  desired  width  2a  of  the  discontinuous  section  to 
the  best  ability  of  the  computational  grid.  The  alternative  of  anchoring  the  potentials  at 
the  two  ends  of  the  well  bottom  would  lead  to  an  approximation  to  a  true  square  potential 
of  width  2a  -  Ax,  with  a  different  corresponding  analytic  transmission  function  for  each 
grid  spacing.  The  choice  of  this  configuration  would  lead  to  unnecessary  error  in  calcula¬ 
tions  at  larger  grid  spacings  Ax,  as  the  well  widths  diverged  from  the  reference  width.  This 
type  of  error  is  visible  in  a  plot  of  transmission  functions  (Figure  5.6),  where  the  period 
of  oscillation  of  the  resonances  in  the  transmission  function  can  be  seen  to  vary  with  the 
grid  spacing,  corresponding  to  changing  well  widths  (80).  By  contrast,  a  similar  series  of 
calculations  made  with  asymmetric  potentials  in  the  style  of  Figure  5.3  shows  the  period 
of  oscillation  remaining  constant  as  the  propagation  time  step  is  varied  (Figure  5.7). 

Aside  from  the  unavoidable  error  incumbent  in  the  discretization  of  a  discontinuous 
potential,  some  error  in  wavepacket  calculations  is  controlled  by  the  selection  of  the  grid 
spacing  Ax  and  the  time  step  At,  based  on  their  conjugate  relationships  with  momentum 
and  energy  quantities,  respectively (7).  Let  the  maximum  available  energy  in  the  system 
be  called 


-E'max  —  Ejint  +  Tmax  +  (5.22) 

where  Eint  is  internal  energy  (not  present  in  these  one-dimensional  cases),  Tmax  =  f?kl^/2n 
is  the  maximum  translational  energy  represented  on  the  grid,  and  Kiax  the  maximum  po¬ 
tential  energy.  The  maximum  momentum,  femax  =  tt/A®,  is  fixed  via  the  FFT  by  its 
conjugate  relationship  with  the  coordinate  grid.  Therefore  a  coordinate  spacing  Axmax  ex¬ 
ists,  beyond  which  the  momentum  grid  will  not  support  the  evolving  wavepackets  needed 
for  the  channel-packet  calculations.  With  regard  to  the  time  step  At,  similar  conditions 
exist  with  regard  to  Specifically,  the  conjugate  relationship  between  energy  and  time 

requires  that  At  <  h/Emax- 
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Probabilicy  of  Reaction 


Energy  (eV) 

Figure  5.6  Convergence  behavior  at  decreasing  grid  spacings  Aa:  for  the  transmission 
functions  of  a  family  of  symmetric  square  wells  whose  width  depends  on  the 
grid  spacing(80). 
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Figure  5.7  Transmission  functions  of  some  asymmetric  square  wells  resembling  Figure 
5.3,  made  by  the  channel-packet  method  in  the  Schrodinger  picture. 


Variable 

Value 

Description 

Xi 

-1 

Lowest  coordinate  in  well 

X2 

1 

Lowest  coordinate  in  right  asymptotic  channel 

-Vo 

-50 

Potential  energy  in  well 

Vi 

0 

Potential  energy  of  left  asymptotic  channel 

V2 

100 

Potential  energy  of  right  asymptotic  channel 

Table  5.1  Potential  parameters  held  constant  for  all  transmission-function  calculations 
in  this  chapter.  All  quantities  are  in  atomic  units. 

The  potentials  used  in  this  chapter  all  share  the  parametric  values  given  in  Table 
5.1.  Grids  are  always  chosen  such  that  the  anchor  points  at  the  leftmost  ends  of  the  well 
and  the  right  asymptote  both  reside  on  the  grid. 

A  baseline  is  needed  to  compare  against  the  results  of  interaction-picture  calculations 
for  which  no  analytic  solutions  exist.  This  role  is  filled  by  a  split-operator  propagator.  The 
split-operator  method  has  become  a  standard  propagation  technique  in  the  Schrddinger  pic¬ 
ture  because  it  is  very  stable  and  delivers  accurate  results  comparatively  quickly(34).  Three 
series  of  calculations  were  run  in  the  Schrodinger  picture  for  comparison  to  interaction- 
picture  results.  All  three  were  transmission  functions  for  various  approximations  to  the 
asymmetric  square  potential  of  Figure  5.2.  The  first  set  fixes  the  coordinate  spacing  at 
Ax  =  0.0025,  and  demonstrates  the  convergence  of  the  transmission  function  toward  the 
analytic  version  with  decreasing  computational  time  step  At.  The  second  set  of  calcula¬ 
tions  uses  a  fixed  time  step  At  =  1.0  •  10”^  atomic  units,  and  demonstrates  convergence 
as  the  coordinate  spacing  Ax  is  decreased.  It  is  recognized  that  changes  in  accuracy  in 
this  series  are  related  to  both  the  decrease  in  the  grid  spacing  and  the  improvement  in  the 
accuracy  of  the  potential  itself  as  finer  grid  spacing  allows  closer  approach  to  truly  vertical 
well  walls.  The  third  set  of  calculations  explores  this  source  of  error  in  two  ways.  First, 
both  At  and  Ax  are  fixed  at  1.0  •  10~®  atomic  units  and  0.01  atomic  units  respectively, 
and  potentials  with  various  degrees  of  slope  in  the  sides  of  the  well  instead  of  the  best 
possible  square-well  approximation  are  tested.  Second,  a  trapezoidal  potential  is  tested 
with  a  fixed  distance  W  =  0.08  atomic  units  between  the  endpoints  of  the  well  bottom 
and  the  beginning  points  of  the  asymptotes,  using  At  =  1.0  •  10“®  atomic  units,  on  grids 
of  various  spacings  Ax  that  represent  the  potential  identically  except  for  the  number  of 
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Variable 

Value 

Description 

“^max 

-10.24 

Lowest  value  on  coordinate  grid  for  Mpller  states 

(J 

0.25 

Width  parameter  for  initial  Gaussians  tpin/out  i^) 

Xq 

0 

Position  parameter  for  “(pin/out  (^) 

ko 

12 

Momentum  parameter  for  t/’m/ouf  (^) 

T 

1 

±0.45 

Mass  parameter  for  -tpin/emt  (^) 

Asymptotic  time  '' 

Table  5.2  Grid  and  wavepacket  parameters  held  constant  for  all  transmission-function 
calculations  in  this  chapter.  All  quantities  are  in  atomic  units. 

points  sampled.  This  series  isolates  the  source  of  degradation  in  accuracy  with  increased 
grid  spacing  to  only  the  grid  itself,  holding  the  potential  constant.  The  convergence  of 
the  transmission  fimction  for  this  series  is  measured  relative  to  the  calculation  with  the 
smallest  Ax.  To  the  greatest  degree  possible,  the  same  parameters  were  used  across  all 
correlation-function  calculations,  in  order  to  isolate  the  variation  in  the  results  to  the  por¬ 
tion  of  the  S-matrix  calculation  that  uses  the  interaction  picture;  namely,  the  derivation 
of  the  Mpller  states. 

The  parameters  listed  in  Table  5.2  are  common  to  all  the  calculations  of  transmission 
function,  in  both  the  SchrOdinger  and  interaction  pictures.  The  grid  parameter  N,  the  total 
number  of  grid  points,  is  varied  in  tandem  with  the  grid  spacing  Ax  to  hold  the  parameter 
a^max  =  AT  Ax  constant  for  all  calculations.  This  ensures  that  sufficient  space  will  exist  on 
the  coordinate  grid  for  the  intermediate  states  '4’in/cmt 

All  calculations  use  a  modified  version  of  Calfas’  code  for  calculation  of  the  correlation 
function,  with  the  parameters  given  in  Table  5.3  common  to  all(80).  The  code  was  modified 
for  this  project  to  support  asymmetric  potentials  and  to  run  on  larger  grids.  The  grid 
spacing  Axe  for  th®  correlation  function  was  always  the  same  as  that  for  the  Mpller 
states,  and  the  number  of  grid  points  Nc  for  the  correlation  function  was  always  twice  the 
number  used  to  calculate  the  Mpller  states.  The  grid  for  the  correlation  function  can  be 
this  small  because  the  code  uses  absorbing  boundary  conditions  of  the  form, 
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Variable 

Value 

Description 

^maxe 

A 

Ate 

ri 

T2 

-20.48 
1.0  •  10-"* 
1.0  •  10"^ 
-0.5 
1.0 

Lowest  value  on  coordinate  grid  for  correlation  function 
Absorbing  boundary  condition  multiplier  (Equation  5.23) 
Time  step  for  correlation  function 

Negative  asymptotic  time 

Positive  asymptotic  time 

Table  5.3  Grid  and  wavepacket  parameters  held  constant  for  all  correlation-function 
calculations  in  this  chapter.  All  quantities  are  in  atomic  units. 


In  this  chapter,  all  calculations  of  the  correlation  function  used  the  value  A  —  1.0  •  10 
given  in  Table  5.3,  along  with  values  of  B  chosen  to  make  \V'  (sTmaxc)!  ~  Such  boundary 
conditions  were  found  for  each  grid  by  trial  and  error  to  prevent  both  reflection  from  and 
transmission  across  the  grid  boundary  by  wavepackets  to  any  detectable  degree.  The  value 
of  Xb  varies  with  the  grid  size.- 

The  calculation  of  the  transmission  function  from  the  correlation  function  can  only 
be  valid  over  a  certain  range  of  energies.  It  can  be  seen  from  equation  2.28  that  this 
method  can  be  numerically  stable  only  over  the  range  of  energies  where  the  wavefunction 
product  rf_  (±A:y)  (±^)  is  numerically  appreciable.  The  energy  spectra  rf^  {E)  of  the 

wavepackets  and  V’out  iised  in  this  chapter’s  calculations  are  shown  with  their  product 
in  Figure  5.8.  Error  in  the  transmission  function  is  measured  as  the  average, 


=  -E 

n  4-rf 


j!CPM  _  ^analytic 


Z=1 


f 


analytic 


(5.24) 


where  the  channel-packet  result  is  compared  to  the  analytic  function  only  at  those 

contiguous  energy  values  where  the  divisor  r]'*  {Ei)  (Ei)  >  0.01.  For  the  wavepackets 
used  here,  this  energy  range  in  atomic  units  is  67  <  E  <  157 . 


5. 2. 3.1  Schrodinger- Picture  Calculations. 


Asymmetric  Square  Wells  With  Various  Time  Steps.  For  the  scenario 
of  this  calculation,  with  a  fixed  Ax  =  0.0025  atomic  units,  the  predicted  maximum  safe 
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Figure  5.8  The  product  r)*  {Ein)  {Eaat)  and  its  components,  the  energy  spectra  of 
the  reactant  and  product  wavepackets.  Since  this  product  is  the  divisor  of 
the  formula  for  the  S  matrix,  its  spectrum  bounds  the  range  of  energies  over 
which  the  S  matrix  may  be  calculated  by  the  channel-packet  method  using  a 
given  pair  of  initial  states  ipin/mit- 
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time  step  is 


At 


max  — 


h 

^max 

h 

Tmax  “1“  Knax 

_ h 

+  Knax 

_ h _ 

(Ax)^  + 

^^ax 

_ 1 _ 

7r2/2  (0.0025)  V  100 

1  •  10"® 


(5.25) 


in  atomic  units.  This  turns  out  to  be  unnecessarily  conservative,  based  as  it  is  on  the 
maximum  possible  total  energy  in  the  model,  rather  than  the  maximum  energy  actually 
seen  in  the  collision.  As  seen  in  Figures  5.9  and  5.10,  accurate  results  are  achieved  over 
the  selected  energy  range  at  much  larger  time  steps.  The  momentum  grid  in  a  square-well 
model  must  necessarily  have  excess  capacity,  since  the  coordinate  spacing  Ax  =  Tr/femax 
must  be  made  small  in  order  to  approximate  the  potential  with  accuracy. 

Figure  5.10  shows  an  order  of  convergence  of  approximately  1.5  at  the  right  end  of 
the  curve,  and  ceases  to  converge  for  time  steps  smaller  than  0.0001  atomic  units.  The 
failure  to  continue  to  converge  results  from  the  inaccuracy  in  the  potential  inherent  to  the 
discrete  representation  with  Ax  =  0.0025  atomic  units. 

Asymmetric  Square  Wells  With  Various  Coordinate  Steps.  A  quadratic 
order  of  convergence  is  evident  when  the  time  step  is  held  constant  and  the  grid  spacing  is 
varied,  as  shown  in  Figure  5.11.  In  this  case  the  convergence  is  limited  by  the  finite  time 
step  At  =  1.0  •  10"®  used  throughout  the  calculations. 

Asymmetric  Trapezoidal  Wells  With  Various  Coordinate  Steps.  The 
previous  result  (Figure  5.11)  includes  two  simultaneously  varying  quantities  that  contribute 
to  error  in  the  transmission  function;  one  being  the  coordinate  spacing  itself,  the  other  the 
variation  of  the  potential  that  is  incumbent  in  seeking  the  best  approximation  to  the 
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Figure  5.9  Transmission  coefficient  of  the  asymmetric  square  well  for  several  choices  of 
propagation  time  step  At  for  the  Mpller  states  only.  The  correlation-function 
step  of  the  channel-packet  process  uses  Ate  =  1-0  *  10  ^  atomic  units  in  all 


cases. 


Bi 


Error 


Figure  5.10  The  error  in  the  transmission  coefficients  of  the  asymmetric  square  well 
transmission  function  as  computed  in  the  Schrodinger  picture.  Error  is 
computed  using  equation  (5.24),  with  the  analytic  square- well  potential  as 
the  reference. 


% 


Figure  5.11  Convergence  of  the  transmission  coefficient  for  the  asymmetric  square  well 
as  the  coordinate  spacing  Ax  is  decreased  in  the  SchrSdinger  picture.  The 
accuracy  of  the  discrete  representation  of  the  square  potential  is  also  im¬ 
proving  with  decreasing  Ax.  Error  is  computed  using  equation  (5.24),  with 
the  analytic  square-well  potential  as  the  reference. 
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Grid  Spacing  of 
“Square”  Well 

Extra  Points  in 
Trapezoidal  Well  Sides 

Equivalent  Ax 

0 

0.01 

1 

0.02 

3 

0.04 

0.08 

7 

0.08 

0.16 

15 

0.16 

Table  5.4  Summary  of  the  potentials  compared  in  Figures  5.12  and  5.17.  The  first  and 
third  columns  are  expressed  in  atomic  units  of  distance.  “Equivalent  Ax"  is 
seen  to  be  idendical  to  the  grid  spacing  of  the  “square"  well  that  is  being 
compared  to  the  trapezoidal  well  on  the  grid  with  spacing  Ax  =  0.01  atomic 
units.  The  trapezoidal  wells  on  the  finer  grid  are  formed  by  adding  enough 
extra  points  to  the  well  sides  (between  the  outer  points  of  the  well  bottom  and 
the  iimer  points  of  the  asymptotes)  that  the  slopes  of  the  well  sides  become 
identical  to  the  corresponding  slopes  of  the  “square"  potentials  on  the  coarser 
grids. 

square  potential  that  is  possible  in  each  grid.  These  factors  can  be  examined  one  at  a 
time  by  mimicking  the  trapezoidal  configuration  of  the  wells  constructed  on  the  more 
coarsely-spaced  grids  on  grids  with  smaller  Ax.  Figure  5.12  holds  the  first  factor  constant, 
comparing  the  convergence  of  trapezoidal  wells,  all  with  At  =  1.0  •  10"®  and  Ax  =  0.01 
atomic  units,  with  various  numbers  of  extra  points  added  to  the  well  walls  to  make  the 
potentials  the  same  shape  as  the  “square”  potentials  previously  constructed  on  grids  with 
larger  values  of  Ax.  For  example,  a  well  with  one  extra  point  between  the  inner  ends  of 
the  asymptotes  and  the  outer  ends  of  the  well  is  compared  to  one  with  no  extra  points 
(a  “square”  well)  on  a  grid  with  Ax  —  0.02  atomic  units.  Table  5.4  enumerates  the 
comparisons  depicted  in  Figure  5.12.  The  figure  shows  that  the  order  of  convergence  to 
the  analytic  square-well  solution  remains  quadratic  when  the  same  potentials  are  modeled 
on  more  finely  spaced  grids,  with  the  finer  grids  showing  somewhat  improved  accuracy  in 
the  transmission  function  relative  to  those  computed  on  the  same  potential  using  a  coarser 
grid. 


Figure  5.13  addresses  the  second  factor.  Here  a  single  trapezoidal  shape  is  maintained 
with  distance  W  =  0.08  atomic  units  between  the  rightmost  point  of  the  well  bottom  and 
the  leftmost  point  of  the  right  asymptote,  while  the  coordinate  spacing  is  varied  over 
the  range  0.0025  <  Ax  <  0.08  atomic  units.  Quadratic  convergence  of  the  transmission 
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Figure  5.12  Error  in  transmission  coefficients  for  asymmetric  trapezoidal  well  potentals 
(dotted  curve),  compared  to  those  for  the  corresponding  asymmtetric  square 
wells  on  coarser  numerical  grids  (solid  curve).  The  trapezoidal  wells  have  no 
grid  points  within  the  well  walls  for  Ax  =  0.01,  one  for  Ax  =  0.02,  three  for 
Ax  =  0.04,  and  so  on,  as  enumerated  in  Table  5.4.  The  entire  channel-packet 
calculation  is  done  in  the  Schrodinger  picture.  Error  is  computed  using 
equation  (5.24),  with  the  analytic  square- well  potential  as  the  reference. 
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Figure  5.13  Convergence  in  the  Schrodinger  picture  of  the  transmission  function  for  a 
trapezoidal  well  of  fixed  shape  as  coordinate  spacing  Ax  is  varied.  The  refer¬ 
ence  transmission  function  uses  Ax  =  0.0025  atomic  units.  All  calculations 
use  At  =  1.0  •  10“®  atomic  units. 

function  is  seen,  not,  of  course,  toward  the  analytic  square- well  transmission  function, 
but  toward  the  most  accurate  calculation  of  the  transmission  function  for  this  particular 
trapezoidal  well;  namely,  the  calculation  performed  using  the  smallest  grid  spacing. 

5.2. 3.2  Interaction- Picture  Calculations. 

Asymmetric  Square  Wells  With  Various  Time  Steps.  For  larger  time 
steps  in  the  Mpller-state  calculation,  the  interaction-picture  approach  shows  linear  conver¬ 
gence  behavior,  at  a  higher  level  of  error  than  in  the  SchrSdinger  picture  for  a  given  time 
step  At.  This  is  illustrated  in  Figure  5.14.  The  interaction-picture  calculation  of  the  Mpller 
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Figure  5.14  Convergence  with  decreasing  time  step  At  of  transmission  functions  com¬ 
puted  in  the  Schrddinger  picture  using  Mpller  states  produced  in  the  in¬ 
teraction  picture,  compared  to  transmission  functions  produced  entirely  in 
the  Schrddinger  picture.  Error  is  computed  using  equation  (5.24),  with  the 
analytic  square-well  potential  as  the  reference. 

states,  while  still  limited  by  the  common  use  of  the  fixed  time  step  Ate  =  1.0  •  10“®  in  the 
Schrodinger  pictiure  for  the  calculation  of  the  correlation  function,  appears  ultimately  to 
converge  to  a  result  of  equal  accuracy  to  that  achieved  by  calculation  of  the  Mpller  states 
in  the  SchrOdinger  picture  at  the  shortest  time  steps.  This  is  reasonably  seen  as  an  effect 
of  the  time  dependence  of  the  interaction-picture  Hamiltonian.  Discontinuous  or  rapidly 
varying  potentials  such  as  square  and  highly  sloped  trapezoidal  wells  and  barriers  would 
be  expected  to  cause  the  interaction-picture  Hamiltonian  to  vary  rapidly  in  time  while  the 
wavepackets  are  in  the  interaction  region  of  the  potential. 
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Asymmetric  Square  Wells  With  Various  Coordinate  Steps.  The  con¬ 
vergence  behavior  of  the  transmission  function  based  on  interaction-picture  Mpller  states  as 
coordinate  spacing  Ax  decreases  is  consistent  with  this  interpretation.  At  the  chosen  small 
value  At  =  Ate  =  1.0  •  10“^  atomic  imits,  Figure  5.14  predicts  equal  accuracy  between 
the  Schrbdinger  and  interaction  pictures.  Figure  5.15  confirms  this  prediction,  with  an 
indication  at  the  largest  grid  spacing,  where  its  Hamiltonian  is  the  least  time-dependent, 
of  somewhat  more  accurate  results  than  the  Schrbdinger  picture  at  the  common  time  step. 
A  similar  experiment  using  the  common  time  step  At  =  Ate  =  l-O  •  10“^  makes  the  point 
even  more  nicely,  as  shown  in  Figure  5.16.  Use  of  the  longer  time  step  exposes  the  inter¬ 
action  picture  to  increased  error  relative  to  the  Schrodinger  picture  for  the  more  steeply 
sloped  potentials,  while  remaining  of  comparable  accuracy  for  the  less  steeply  sloped  po¬ 
tentials.  Constraints  on  the  coordinate  grid  required  to  keep  the  potential’s  anchor  points 
at  X  =  ±1  atomic  imit  on  the  grid  restrict  the  available  number  of  large  values  of  Ax  be¬ 
low  the  momentum-grid  limit.  Both  the  interaction- picture  and  the  Schrodinger  methods 
appear  to  converge  toward  the  analytic  result  approximately  quadratically  in  this  scenario. 

Asymmetric  Trapezoidal  Wells  With  Various  Coordinate  Steps.  The 
interaction  picture  does  not  show  the  same  benefit  from  decreasing  the  coordinate  spacing 
as  the  Schrodinger  picture  in  the  test  used  here,  where  the  time  step  is  held  constant  and 
the  coordinate  spacing  is  held  constant  at  Ax  =  0.01  atomic  units,  while  the  slope  of  the 
well  sides  is  varied.  The  interaction  picture  remains  similar  in  accuracy  to  the  SchrSdinger 
picture  at  this  time-step  size,  as  shown  in  Figure  5.17. 

Figure  5.18  demonstrates  that  the  interaction  picture  converges  similarly  toward  the 
benchmark  Schrodinger-picture  calculation  of  the  transmission  function  for  the  fixed-shape 
trapezoidal  well  as  grid  spacing  decreases.  Not  surprisingly,  the  approach  of  the  interaction- 
picture  version  to  the  benchmark  ceases  at  higher  grid  resolutions  as  the  interaction  picture 
converges  to  its  own  best  version  of  the  transmission  function. 
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Error  with  varying  grid  spacing  Ax  of  transmission  functions  for  the  asym- 
metric  square  well  with  M0ller  states  computed  in  the  interaction  picture, 
compared  to  results  obtained  entirely  in  the  SchrOdinger  picture.  These  cal¬ 
culations  use  At  =  Ate  =  1-0  •  10'^  atomic  units.  Error  is  computed  using 
equation  (5.24),  with  the  analytic  square- well  potential  as  the  reference. 


0.001 


0.01 


A  X  {a.u} 


Figure  5.16  Error  with  varying  grid  spacing  Arc  of  transmission  functions  for  the  asym¬ 
metric  square  well  with  Mpller  states  computed  in  the  interaction  picture, 
compared  to  results  obtained  entirely  in  the  SchrSdinger  picture.  These  cal¬ 
culations  use  At  =  Ate  =  1.0  •  10"^  atomic  units.  Error  is  computed  using 
equation  (5.24),  with  the  analytic  square- well  potential  as  the  reference. 
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Error 


Figure  5.17  Frror  in  transmission  functions  for  asymmetric  trapezoidal  well  potentals 
(dotted  curve),  compared  to  those  for  the  corresponding  asymmtetric  square 
wells  on  coarser  munericaJ  grids  (solid  curve).  The  trapezoidal  wells  have 
no  grid  points  within  the  well  walls  for  Ax  =  0.01,  one  for  Ax  =  0.02,  three 
for  Ax  =  0.04,  and  so  on,  as  enumerated  in  Figure  5.4.  The  Mpller  states 
for  these  channel-packet  calculations  are  obtained  in  the  interaction  pic¬ 
ture.  The  equivalent  curve  for  trapezoidal  wells  modeled  in  the  Schrodinger 
picture  (dashed  line)  is  included  for  comparison.  Error  is  computed  using 
equation  (5.24),  with  the  analytic  square- well  potential  as  the  reference. 
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Figure  5.18  Convergence  in  the  interaction  picture  of  the  transmission  function  for  a 
trapezoidal  well  of  fixed  shape  as  coordinate  spacing  Ax  is  varied.  The  refer¬ 
ence  transmission  function  uses  Ax  =  0.0025  atomic  units  in  the  SchrSdinger 
picture.  All  calculations  use  At  =  1.0  •  10'^  atomic  units. 


5.3  Summary 

The  SchrOdinger-picture  split-operator  propagator  delivers  results  that  improve  in 
accuracy  consistently  with  increasing  coordinate  and  time  sampling  rates.  It  can  be  used 
as  a  basis  for  comparison  for  the  interaction-picture  propagator.  The  calculation  of  Mpller 
states  in  the  interaction  picture  is  also  a  valid  technique,  but  must  take  into  account  the 
rapidity  of  changes  in  the  potential  in  the  choice  of  time  steps  in  order  to  achieve  accuracy 
similar  to  that  attained  by  the  Schrodinger-picture  method. 
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VI.  Application  of  the  Interaction  Picture  in  Reactive  Scattering 

This  chapter  uses  Mpller-state  and  S- matrix  calculations  on  various  one-dimensional  poten¬ 
tials  to  demonstrate  the  utility  of  the  interaction  picture  in  reactive  scattering  calculations. 
Various  potentials  of  the  general  form  (Figure  6.1), 

V  (x)  =  -t-  -I-  DQ  (x  —  a)  0  (c  —  x)  x  -I-  VqQ  (x  —  a) , 

(6.1) 

where  0  (x)  represents  a  Heaviside  step  function,  are  chosen  to  schematically  represent 
the  reaction  path  of  a  reactive  molecular  collision  that  may  have  two  different  asymptotic 
Hamiltonians.  A  potential  barrier  of  the  form 

V(x)  =  Asech(Bx)  (6.2) 

is  used  to  demonstrate  a  type  of  case  where  the  interaction  picture  is  advantageous  relative 
to  the  Schr5dinger  pictmre  for  the  channel-packet  technique. 

6.1  Application  to  a  Reactive  Potential 

The  channel-packet  method’s  applicability  to  potentials  with  multiple  asymptotic 
potential  energies  was  demonstrated  in  both  the  interaction  and  Schrodinger  pictures  in 
Chapter  V.  For  this  potential,  as  x  ^  — oo,  the  asymptotic  Hamiltonian  is 


Hi  =  pV2/x, 


(6.3) 


and  as  X  — >  — cx),  the  asymptotic  Hamiltonian  is 

Hf  =  pV2M  +  V^).  (6.4) 

An  example  using  the  values  for  the  parameters  in  equation  (6.1)  given  in  Table  6.1  appears 
as  the  dotted  curve  in  Figure  6.2. 
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V  ( X )  { atomic  units } 


Figure  6.1  An  asymmetric  triple  Gaussian  potential. 


Peurameter 

Value 

A 

0.03 

B 

0.09 

C 

0.02 

D 

0.00125 

a 

-4.0 

b 

0.0 

c 

4.0 

a 

1.0 

1.0 

7 

1.0 

Table  6.1  The  coefficients  used  in  Equation  6.1  to  create  the  potential  function  that 
appears  in  Figures  6.2  et.  seq.  All  quantities  use  atomic  units. 


Parameter 

Value 

Xq 

0.0 

feo 

8.5 

G 

0.55 

1224 

Table  6.2  The  coefficients  used  in  Equation  5.3  to  create  the  asymptotic  wavepackets 
shown  in  Figure  6.2.  All  quantities  use  atomic  imits. 

The  transmission  function  of  this  potential  can  be  computed  by  the  channel-packet 
method  using  the  interaction  picture,  beginning  with  the  product  and  reactant  wavepack¬ 
ets  shown  in  Figure  6.2.  The  choice  of  -|-A:o  in  equation  (5.3)  will  yield  the  probability,  as 
a  function  of  kinetic  energy,  that  reactants  approaching  from  the  left  will  form  products 
exiting  to  the  right.  The  values  of  the  coefficients  used  in  equation  (5.3)  used  to  generate 
the  reactant  and  product  states  are  given  in  Table  6.2.  In  the  SchrOdinger  picture,  the  first 
propagation  in  equation  (2.18)  is  performed  analytically (27),  using  the  initial  wavepack¬ 
ets  together  with  and  to  obtain  intermediate  reactant  and  product  wavepackets 
respectively.  As  shown  in  Figure  6.3,  these  wavepackets  undergo  both  translation  and 
spreading  relative  to  the  initial  wavepackets. 

Typically,  this  wavepacket  translation  and  spreading  in  the  SchrSdinger  picture  gen¬ 
erates  a  requirement  for  large  grids.  However,  in  the  interaction  picture,  the  intermediate 
wavepackets  do  not  translate  or  spread,  and  remain  identical  to  the  initial  wavepackets. 
Figure  6.3  also  illustrates  the  intermediate  interaction-picture  chaimel  packets,  and  demon- 
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Probability  Density 


Initial  Channel  Packets 


t  =  0  atomic  units 


X  { atomic  units } 


Figure  6.2  The  initial  channel  packets  (x)  and  •0ou4  (^)  (solid  line),  are  the  same  in 
the  interaction  and  Schrddinger  pictures,  since  they  are  evaluated  at  t  =  0. 
The  potential  (dotted  line)  is  the  sum  of  two  Gaussian  barriers,  a  Gaussian 
well,  and  a  ramp  function  which  is  zero  for  x  <  —4,  0.01  fora:  >  4,  and  rises 
linearly  in  between. 
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V  ( X )  { atomic  units } 


Intermediate  States 
t=  ±2000  atomic  units 


Figure  6.3  The  reactant  channel  packet  propagated  backward  in  time  to  t  =  —2000 
atomic  units,  and  the  product  channel  packet  propagated  forward  in  time  to 
t  =  2000  atomic  units.  Since  the  propagation  occurs  under  a  free-particle 
Hamiltonian,  the  wavepackets  are  unaffected  in  the  interaction  picture  (solid 
line).  The  Schrbdinger-picture  packet  (dashed  lines)  translates  to  the  left  (re¬ 
actant,  long  dashes)  or  right  (product,  short  dashes),  and  spreads,  requiring 
a  larger  grid. 
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Parameter 

Value 

T 

±2000 

to 

8.5 

At  I 

0.5 

Ats 

1.0 

Table  6.3  The  quantities  associated  with  the  generation  of  the  intermediate  states  in 
Figure  6.3  and  the  Mqller  states  in  Figure  6.4.  Intermediate  states  (Figure 
6.3)  are  calculated  analytically  for  the  Gaussian  asymptotic  states  at  time 
T  =  —2000  atomic  units  for  the  reactant  state,  and  r  =  +2000  atomic  units 
for  the  reactant  state.  Mpller  states  (Figure  6.4)  are  generated  by  propagating 
the  intermediate  states  to  the  time  to  =  0,  using  the  computational  time  step 
At  I  =  0.5  atomic  units  in  the  interaction  picture,  and  Ats  =  1.0  atomic  units 
in  the  SchrOdinger  picture. 

strates  that  they  require  a  smaller  grid  when  compared  with  the  intermediate  Schrodinger- 
picture  wavepackets.  In  the  interaction  picture,  the  intermediate  wavepackets  are  used 
in  equation  (3.36)  together  with  the  appropriate  asymptotic  Hamiltonian  to  compute  the 
reactant  and  product  Mpller  states.  The  Mpller  states  shown  in  Figure  6.4  were  computed 
on  a  grid  of  256  points  using  the  nested  interaction  picture  with  a  four-dimensional  Krylov 
subspace.  Constants  used  for  the  propagation  are  listed  in  Table  6.3.  For  comparison,  the 
same  Mpller  states  are  computed  in  the  Schrbdinger  picture,  requiring  a  grid  of  512  points. 

To  complete  the  calculation,  the  Mpller  states  were  propagated  in  the  Schrodinger 
picture  using  a  spht-operator  propagator  together  with  absorbing  boundary  conditions 
to  compute  the  correlation  function  in  equation  (2.26),  represented  in  Figure  6.5(23,80). 
S-matrix  elements  are  then  computed  using  equation  (2.28),  resulting  in  the  probability 
for  reaction  shown  in  Figure  6.6.  For  comparison,  the  probability  of  reaction  computed 
entirely  within  the  Schrodinger  picture  is  also  shown  in  Figure  6.6.  It  is  important  to 
note  that  since  short  iterative  Lanczos  propagation  is  employed  when  using  the  nested 
interaction  picture,  a  greater  number  of  FFTs  is  required  per  time  step  when  compared 
to  the  split-operator  approach  commonly  used  for  time-independent  Hamiltonians  in  the 
Schrodinger  picture. 


6-6 


Probability  Density 


Figure  6.4  The  reactant  Moller  state,  (solid  line),  is  the  result  of  propagating  the 

intermediate  reactant  state  forward  in  time  to  t  =  0,  where  the  interaction 
and  SchrSdinger  pictures  are  again  identical.  The  dashed  line  is  the  product 
Mpller  state,  the  result  of  propagating  the  intermediate  product  state 

backward  in  time  to  t  =  0. 
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Probability  Density 


Correlation  Function 


Time  { atomic  units } 


Figure  6.5  The  absolute  value  squared  of  the  correlation  function,  computed  using  equa¬ 
tion  (2.26). 


Probability  of  Transmission 


Parameter 

Potential  1 

Potential  2 

Potential  3 

Potential  4 

A 

0.03 

0.03 

0.03 

0.03 

B 

0.09 

0.09 

0.09 

0.09 

C 

0.03 

0.03 

0.03 

0.03 

D 

0.0 

0.0 

0.0 

0.0 

a 

-4.0 

-4.0 

-4.0 

-21.0 

b 

0.0 

0.0 

0.0 

0.0 

c 

4.0 

4.0 

4.0 

21.0 

a 

1.0 

0.5 

0.25 

0.015625 

0 

1.0 

0.5 

0.25 

0.015625 

7 

1.0 

0.5 

0.25 

0.015625 

Table  6.4  The  coefficients  used  in  Equation  6.1  to  create  four  potential  functions  with 
similar  energy  characteristics  but  differing  slopes.  All  quantities  use  atomic 
units. 

6.2  Effect  of  Potential  Slope 

The  interaction  picture  is  at  a  performance  disadvantage  relative  to  the  SchrOdinger 
picture  regarding  time-independent  potentials.  Since  all  spatially  non-constant  potentials 
yield  time-dependent  Hamiltonians  in  the  interaction  picture,  the  length  of  the  calcula¬ 
tion  time  step  At  is  more  constrained  than  it  is  for  the  corresponding  time-independent 
Hamiltonian  in  the  SchrOdinger  picture.  In  this  section,  the  effect  of  varying  the  rapidity 
of  change  of  the  potential  with  the  spatial  coordinate  is  examined,  using  several  symmetric 
potentials  (Figure  6.7)  fomulated  according  to  equation  (6.1),  with  the  parameters  given  in 
Table  6.4.  The  measure  of  error  chosen  is  the  amplitude  error  of  the  wave  function,  given 
by  equation  4.2.  For  a  fixed  time  step,  chosen  to  be  At  =  1.0  atomic  units,  three  reactant 
Mpller  state  calculations  are  performed  on  identical  coordinate  grids.  The  propagation 
techniques  examined  are  a  Lanczos  method  in  the  interaction  picture  using  a  first-order 
approximation  of  the  Hamiltonian,  a  second-order  Lanczos  method  in  the  interaction  pic¬ 
ture,  and  a  split-operator  method  in  the  Schrodinger  picture.  The  amplitude  error  of  each 
is  measured  periodically  relative  to  a  Schrodinger-picture  split-operator  propagation  using 
a  time  step  At'  =  0.1.  Before  comparison  to  the  reference  wavepacket,  the  interaction- 
picture  wavepackets  of  course  must  be  converted  to  the  SchrOdinger  picture. 

The  calculation  on  the  steepest  potential  (Potential  1)  shows  that  for  this  poten¬ 
tial,  the  chosen  time  step  produces  accurate  Mpller  states  from  both  the  Schrodinger  and 
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Figure  6.7  The  four  triple  Gaussian  potentials  used  to  investigate  the  effect  of  varying 
potential  slope  on  the  accuracy  of  propagation  in  the  interaction  picture. 
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Figure  6.8  The  growth  in  amplitude  error  (equation  (4.2))  during  the  computation  of 
the  reactant  Mpller  state  for  Potential  1.  The  solid  line  corresponds  to  a 
scond-order  nested  interaction-picture  Lanczos  propagator,  the  alternating 
dots  and  dashes  to  a  first-order  one.  The  dotted  line  shows  the  error  in  the 
computation  of  the  same  Mpller  state  using  a  split-operator  propagator  in 
the  Schrodinger  picture.  All  propagations  shown  used  a  time  step  At  =  1.0 
atomic  unit.  The  amplitude  error  for  all  three  propagations  is  measured 
against  a  split-operator  propagation  using  At'  =  0.1  atomic  units. 
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the  second-order  interaction  pictures,  as  shown  in  Figure  6.8.  The  first-order  interaction- 
picture  propagator  diverges  rapidly  and  fails  to  produce  an  acceptable  Mpller  state.  The 
second-order  interaction-picture  propagator  diverges  as  the  wavepacket  encounters  the  po¬ 
tential,  but  maintains  good  correspondence  to  the  reference  calculation  overall.  The  much 
more  gradual  divergence  of  the  Schrodinger-picture  method  begins  to  accelerate  late  in  the 
calculation,  but  still  results  in  a  M0ller  state  that  is  markedly  closer  to  the  reference  than 
the  interaction-picture  state  is.  This  is  reasonable  to  expect,  since  the  reference  state  is 
computed  by  the  identical  method,  the  only  difference  being  the  length  of  the  time  step. 

The  somewhat  reduced  slope  of  Potential  2,  with  its  slightly  broader  features,  worked 
to  the  benefit  of  all  three  propagators,  though  not  enough  to  get  an  accurate  result  from 
the  first-order  interaction-picture  proi)agator.  The  divergence  of  both  the  second-order 
interaction-picture  and  the  Schrbdinger-picture  propagators  is  less  rapid,  and  occurs  later, 
even  though  the  interaction  region  of  the  potential  is  encountered  earlier  in  the  propagation. 
Figure  6.9  depicts  the  evolution  of  the  M0ller-state  error  for  Potential  2. 

Figure  6.10  shows  the  error  associated  with  M0ller-state  propagation  on  Potential  3, 
which  has  still  broader  and  more  gradually  sloped  features  than  Potentials  1  and  2.  The 
marginal  improvement  in  the  results  from  all  three  propagations  seen  between  Potentials 
1  and  2  is  extended  with  Potential  3. 

Potentials  1  through  3  are  all  very  similar,  and  demonstrate  the  expected  benefits 
of  reduced  potential  slope  to  the  interaction  picture,  while  showing  that  the  SchrOdinger 
picture  also  enjoys  improved  accuracy.  To  achieve  yet  broader  and  more  gradually  sloped 
features  with  a  triple  Gaussian  potential,  while  retaining  the  well  and  barrier  energies, 
requires  moving  the  locations  of  the  outer  Gaussians  farther  away  from  the  central  one. 
Potential  4  is  the  result  of  such  an  operation.  The  error  of  the  M0ller-state  calculations 
on  this  potential,  shown  in  Figure  6.11,  is  higher  than  might  be  expected  from  the  trend 
seen  with  Potentials  1  through  3;  however,  the  larger  extent  of  this  potential  necessitated 
four  times  larger  grids  in  both  the  interaction  and  Schrodinger  pictures.  The  problem  is 
therefore  of  somewhat  greater  computational  difficulty. 
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Figure  6.9  The  growth  in  amplitude  error  (equation  (4.2))  during  the  computation  of 
the  reactant  Mpller  state  for  Potential  2.  The  solid  line  corresponds  to  a 
scond-order  nested  interaction-picture  Lanczos  propagator,  the  alternating 
dots  and  dashes  to  a  first-order  one.  The  dotted  line  shows  the  error  in  the 
computation  of  the  same  Mpller  state  using  a  split-operator  propagator  in 
the  SchrOdinger  picture.  All  propagations  shown  used  a  time  step  At  =  1.0 
atomic  imit.  The  amplitude  error  for  all  three  propagations  is  measiured 
against  a  split-operator  propagation  using  At'  —  0.1  atomic  units. 
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Amplitude  Error 


Figure  6.10  The  growth  in  amplitude  error  (equation  (4.2))  during  the  computation  of 
the  reactant  Mpller  state  for  Potential  3.  The  solid  line  corresponds  to  a 
scond-order  nested  interaction-picture  Lanczos  propagator,  the  alternating 
dots  and  dashes  to  a  first-order  one.  The  dotted  line  shows  the  error  in  the 
computation  of  the  same  Mpller  state  using  a  split-operator  propagator  in 
the  Schrddinger  picture.  All  propagations  shown  used  a  time  step  At  =  1.0 
atomic  unit.  The  amplitude  error  for  all  three  propagations  is  measured 
against  a  split-operator  propagation  using  At'  =  0.1  atomic  units. 
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Amplitude  Error 


Time  { atomic  units } 


Figure  6.11  The  growth  in  amplitude  error  (equation  (4.2))  during  the  computation  of 
the  reactant  Mpller  state  for  Potential  4.  The  solid  line  corresponds  to  a 
scond-order  nested  interaction-picture  Lanczos  propagator,  the  alternating 
dots  and  dashes  to  a  first-order  one.  The  dotted  line  shows  the  error  in  the 
computation  of  the  same  Mpller  state  using  a  split-operator  propagator  in 
the  Schrodinger  picture.  All  propagations  shown  used  a  time  step  At  =  1.0 
atomic  imit.  The  amplitude  error  for  all  three  propagations  is  measured 
against  a  split-operator  propagation  using  At'  =  0.1  atomic  units. 
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Figure  6.12  places  the  error  trajectories  for  the  second-order  interaction-picture  prop¬ 
agator  with  all  four  potentials  together.  The  general  conclusion  to  be  drawn  is  that  error 
levels,  while  less  for  lower-sloped  potentials  on  the  same  computational  grid,  are  not  sig¬ 
nificantly  altered  when  the  barrier  and  well  energies  remain  the  same. 


Figure  6.12  The  amplitude  error  functions  for  the  second-order  interaction-picture  prop¬ 
agator  on  the  four  symmetric  triple  Gaussian  potentials,  displayed  together 
for  comparison. 

6.3  Effect  of  Wavepacket  Compactness 

Employment  of  the  interaction  picture  in  the  calculation  of  M0ller  states  can  be  re¬ 
lied  upon  to  enable  reduction  of  the  computational  coordinate  grid  by  a  factor  of  two, 
as  demonstrated  above.  However,  a  factor  of  two  is  not  enough  to  result  in  a  reduction 
in  the  propagation  time  relative  to  the  Schrbdinger  picture.  This  is  because  the  (usually 
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shorter)  time  steps  in  the  Lanczos  method  require  at  least  as  many  FFTs  as  there  are 
Krylov  basis  vectors,  as  compared  to  two  FFTs  per  time  step  in  the  split-operator  propa¬ 
gation  scheme,  which  can  be  employed  with  the  time-independent  Hamiltonians  common 
to  the  Schrodinger  picture.  The  nested  interaction  picture  requires  an  additional  pair  of 
FFTs  per  time  step  to  compute  the  expectation  values  (x)  and  (p) ;  and  the  second-order 
Lanczos  algorithm  adds  an  expensive  matrix  diagonalization  as  well  in  order  to  lengthen 
the  allowable  time  step.  This  is  a  worthwhile  investment,  as  shown  in  Section  6.2  where 
results  of  a  first-order  and  a  second-order  propagator  are  compared. 

The  Schrodinger  picture  can  require  much  larger  grids  than  the  interaction  picture 
when  the  wavepackets  remain  compact  throughout  the  propagation.  A  wavepacket  whose 
momentum  has  a  large  absolute  value  requires  a  large  value  of  fcmax,  with  a  concomitant 
small  coordinate  spacing  Ax.  This  in  turn  necessitates  a  large  number  N  of  grid  points 
to  contain  the  trajectory  of  the  coordinate  representation,  even  though  the  wavepackets 
in  both  representations  may  be  relatively  compact.  However,  compact  wavepackets  my  be 
propagated  in  the  nested  interaction  picture  using  grids  barely  large  enough  to  contain  the 
wavepackets  alone;  the  grids  are  adjusted  continuously  to  follow  their  trajectories  in  both 
representations.  This  situation  is  illustrated  with  the  barrier  potential  shown  in  Figure 
6.13.  Two  collision  scenarios  are  examined  with  this  potential,  both  at  a  kinetic  energy 
of  3.0  atomic  units.  One  scenario  involves  a  reactant  state  with  low  reduced  mass  and 
correspondingly  low  momentum;  the  other  a  reactant  state  with  higher  reduced  mass  and 
momentum. 

6.S.1  Low-Momentum  Collision.  For  this  example,  a  reactant  state  with  reduced 
mass  /i  =  18.36  atomic  units  is  chosen,  with  momentum  atomic  units  and  Gaussian  width 
parameter  a  =  0.05.  The  propagation  begins  in  the  asymptotic  region  at  time  — r  = 
—20  atomic  units,  and  proceeds  in  time  steps  At  =  0.01  atomic  units.  Mpller  states  are 
derived  in  both  the  interaction  and  SchrSdinger  pictures.  The  results  of  this  scenario 
are  summarized  in  Table  6.5.  An  eightfold  reduction  in  the  grid  size  required  for  the 
calculation  was  realized  in  the  interaction  picture  as  compared  to  the  Schrbdinger  pictme. 
This  is  sufficient  for  the  interaction  picture  to  be  the  faster  of  the  two  computational 
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V  (x)  { atomic  units } 


Figure  6.13  The  barrier  potential  V  (x)  =  3.0  sech  (x). 
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Propagation 

Method 

Grid 
Size  (JV) 

Grid 

Spacing  (Ax) 

CPU 

Time  (sec) 

SP/SPO 

512 

0.1 

6.8 

SP/SPO 

256 

3.2 

IP/2SIL 

32 

2.9 

Table  6.5  M0ller-state  propagations  with  a  low-momentum  reactant  state  on  the  bar¬ 
rier  potential  of  Figure  6.13,  using  the  split-operator  (SPO)  method  in  the 
Schrodinger  picture  (SP)  and  a  second-order  short  iterative  Lanczos  (2SIL) 
method  in  the  interaction  picture  (IP)  . 


Propagation 

Method 

Grid 
Size  (N) 

Grid 

Spacing  (Ax) 

CPU 

Time  (sec) 

SP/SPO 

1024 

0.1 

14.8 

SP/SPO 

512 

0.1 

6.8 

IP/2SIL 

32 

0.2 

2.9 

Table  6.6  Mpller-state  propagations  with  a  low-momentum  reactant  state  on  the  bar¬ 
rier  potential  of  Figure  6.13,  using  the  split-operator  (SPO)  method  in  the 
Schrodinger  picture  (SP)  and  a  second-order  short  iterative  Lanczos  (2SIL) 
method  in  the  interaction  picture  (IP)  . 

techniques  in  this  scenario.  The  CPU  times  reported  in  the  table  are  from  the  Silicon 
Graphics  IRIX  “time”  utility  on  MIPS  10000  processors. 

6.3.2  High-Momentum  Collision.  A  higher-momentum  wavepacket  makes  more 
clear  a  weakness  of  the  Schrodinger  picture  when  compared  to  the  interaction  picture.  At 
the  same  kinetic  energy  as  the  previous  case,  if  the  reactant  state  has  mass  fi  =  1836 
atomic  units  and  momentum  ko  =  104.2,  a  starting  time  — r  =  —100  is  required  for 
the  intermediate  wavepacket  to  be  clear  of  the  interaction  region  at  the  beginning  of  the 
forward  propagation.  The  greater  time  than  the  low-momentum  case  results  from  more 
rapid  spreading  of  the  wavepacket  during  its  analytic  propagation  backward  in  time  to 
create  the  intermediate  state  at  t  =  — r.  The  end  result  is  a  need  for  more  room  in  the 
coordinate  representation  to  contain  the  extra  spreading  and  translation  of  the  wavepacket. 
In  contrast,  the  interaction-picture  propagation  can  be  done  on  the  same  grid  as  the 
previous  calculation.  Table  6.6  summarizes  the  results  of  this  group  of  propagations,  which 
use  a  time  step  At  =  0.05  atomic  units  on  MIPS  10000  processors. 
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6-4  Summary 

Propagation  in  the  interaction  picture  is  a  reliable,  accurate  option  for  computing  the 
M0ller  states  needed  for  the  channel-packet  method  of  deriving  selected  S-matrix  elements. 
The  interaction  picture  can  usually  accomplish  the  calculation  on  a  grid  that  is,  at  most, 
half  the  minimum  size  possible  when  propagating  in  the  Schrodinger  picture.  However, 
because  of  the  larger  number  of  FFTs  and  matrix  diagonalizations  per  time  step,  as  well 
as  the  time  dependence  of  its  Hamiltonians,  propagation  in  the  interaction  picture  tends 
to  require  more  time  to  reach  an  accurate  evaluation  of  a  Mpller  state  than  propagation 
in  the  SchrSdinger  picture.  The  interaction  picture  becomes  most  advantageous  in  grid 
savings  and  comparative  computational  time  in  high-momentum  collisions. 
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VII.  Conclusion 


Through  the  use  of  the  channel-packet  method,  the  interaction  picture  can  be  used  to 
compute  reactive  scattering  matrix  elements  on  reduced  computational  grids.  Previous 
use  of  the  interaction  picture  has  been  restricted  to  single  scattering  channels,  limiting 
its  application  to  non-reactive  scattering  problems.  The  channel-packet  method  opens 
the  possibility  of  extending  the  efficacy  of  the  interaction  picture  by  allowing  each  Mpller 
state  to  be  computed  in  an  interaction  picture  derived  from  its  asymptotic  Hamiltonian.  A 
twofold  reduction  in  required  grid  size  is  generally  possible,  and  under  the  right  conditions, 
larger  grid  reductions  can  be  achieved. 

7. 1  Efficacy  of  the  Interaction  Picture  in  One  Dimension 

Three  primary  achievements  have  been  demonstrated  during  this  research  project. 
First,  a  method  has  been  found  to  apply  the  interaction  picture  to  multichannel  reactive 
scattering,  retiring  the  commonplace  that  “there  is  no  convenient  generalization  of  the 
interaction  picture”  to  such  scenarios  (1).  Second,  the  nested  interaction  picture  has 
been  shown  reliably  to  allow  computation  of  Mpller  states  on  grids  at  least  a  factor  of  two 
smaller  than  the  smallest  possible  grid  required  to  derive  the  same  states  in  the  SchrOdinger 
picture.  Third,  while  accurate  Mpller  states  for  most  one-dimensional  potentials  can  be 
computed  more  quickly  in  the  Schrbdinger  picture,  under  conditions  where  a  Mpller  state 
that  remains  compact  in  both  its  coordinate  and  its  momentum  representation  allows  much 
larger  grid  size  reductions,  the  nested  interaction  picture  has  been  shown  to  be  the  faster 
choice  in  one  dimension.  What  happens  to  this  situation  in  higher-order  systems  remains 
to  be  shown,  but  as  suggested  in  Section  7.2,  the  advantages  of  the  Schrbdinger  picture 
may  erode  as  more  degrees  of  freedom  are  added. 

The  interaction  picture  is  now  indisputably  applicable  to  the  calculation  of  reactive  S- 
matrix  elements,  and  has  been  shown  to  dovetail  effectively  with  the  channel-packet  method 
and  absorbing  boundary  conditions  in  reducing  the  memory  required  in  the  computation 
of  Mpller  states.  Because  of  the  time-dependence  of  the  interaction-picture  Hamiltonian, 
however,  the  interaction  picture  requires  shorter  time  steps  in  the  interaction  region  of  the 


Method. 

Comment 

Sequential 

Finite  Basis 
Sequential  Nested 
“Heisenberg”  Nested 

Reliable,  accurate,  slow,  no  grid  reduction 

Unstable,  unsuitable  for  long  propagations 

Reliable,  accurate,  slow,  grid  reduction  at  least  twofold 
Needs  analytic  expression  for  potiential;  works  only  for 
certain  potentials 

Table  7.1  Summary  of  results  of  investigations  of  various  interaction-picture  techniques. 

potential,  so  in  one  dimension,  reductions  in  computational  time  relative  to  the  Schrodinger 
picture  are  realized  only  under  special  circumstances,  with  the  largest  reductions  in  grid 
size.  Of  the  numerous  approaches  to  wavepacket  propagation  in  the  interaction  picture 
investigated  in  this  project,  only  the  laborious  sequential  methods  were  reliably  accurate 
when  implemented  and  tested,  and  only  the  sequential  nested  interaction  picture  was  able 
to  deliver  any  reduction  in  grid  size  requirements.  Table  7.1  summarizes  the  methods  tested 
and  the  results  in  each  case.  All  implementations  used  the  iterated  Lanczos  algorithm  to 
represent  the  time-evolution  operator. 

7.S  Extension  of  the  Nested  Interaction  Picture  to  More  Degrees  of  Freedom 

If  the  interaction  picture  is  ever  to  be  useful  in  genuine  molecular  scattering  calcu¬ 
lations,  it  must  be  applied  in  three  or  more  dimensions.  The  interaction  picture  offers  a 
reduction  in  the  memory  requirements  for  such  calculations  by  reducing  the  grid  require¬ 
ments  along  the  translational  coordinates  of  a  multidimensional  model.  Computational 
time,  a  more  desirable  commodity  to  conserve  than  simply  memory,  is  not  necessarily  re¬ 
duced  in  one-dimensional  calculations  using  the  interaction  picture  because  of  the  shorter 
time  steps  necessitated  by  the  interaction  picture’s  time-dependent  Hamiltonians.  How¬ 
ever,  with  added  dimensions,  the  nature  of  the  underlying  mathematics  suggests  that 
reductions  in  computational  time  could  be  realized  in  multidimensional  calculations  even 
if  the  one-dimensional  interaction  picture  process  requires  more  time  than  the  Schrodinger 
picture  on  a  larger  grid.  Assuming  that  the  FFT  dominates  the  computational  effort  re¬ 
quired  in  the  calculation  of  the  M0ller  states,  effort  En  scales  with  number  n  of  degrees  of 
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freedom  as 


En  =  (n-  l)!P„logP„, 


(A.7) 


where  the  argument  P„  is  the  product  of  the  grid  sizes  for  all  n  degrees  of  freedom.  The 
derivation  of  equation  (A.7)  appears  in  Appendix  A. 

Suppose  the  nested  interaction  picture  in  n  dimensions  allows  one  dimension’s  coor¬ 
dinate  grid  to  be  reduced  from  Nj  to  Nj/b.  If  the  interaction-picture  propagation  requires 
a  times  longer  than  the  Schrodinger  picture  in  one  dimension,  then  using  equation  (A.7), 


a  = 


fTA 

aNjlogNj 

ablogNj 

£(i.  £iL\ 

ab\  logNjJ' 


(7.1) 


for  some  pair  of  real  constants  a  and  /?,  which  account  for  differences  in  the  optimum  time 
step  and  in  the  number  of  FFTs  needed  per  time  step  between  the  two  methods.  The  ratio 
of  propagation  times  for  n  >  1  dimensions  is 


\TsJn  aP„logP„ 

1  logfe  {P-a)^\og^ 

b  61ogP„  aP„logPn 


As  the  number  of  degrees  of  freedom  n  increases,  the  terms  with  log  Pn  in  the  denominator 
vanish,  and  the  computational-effort  ratio  given  by  equation  (7.2)  approaches 


lim 

n— >oo 


1 

b' 


(7.3) 


Hence,  as  degrees  of  freedom  are  added,  the  comparative  computational  effort  of  the  in¬ 
teraction  picture  relative  to  the  Schrbdinger  picture  approaches  the  ratio  of  the  required 
grid  sizes  in  the  single  dimension  that  is  different. 
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1.3  Topics  for  Further  Study 

The  successes  obtained  within  the  confines  of  this  investigation  point  the  way  to 
many  related  topics  where  the  use  of  interaction  picture  continues  to  be  an  attractive  po¬ 
tential  tool  for  scattering  calculations.  The  unsatisfactory  performance  of  the  finite-basis 
and  “Heisenberg”  approaches  has  not  been  conclusively  explored.  Either  of  these  tech¬ 
niques  would  be  a  significant  advance  over  what  has  thus  far  been  demonstrated  should 
the  problems  be  overcome.  The  finite-basis  approach  may  well  be  inherently  unstable  be¬ 
cause  it  bases  its  Krylov  vectors  on  a  representation  of  the  Hamiltonian  in  another  Krylov 
subspace,  instead  of  in  either  the  coordinate  or  momentum  representation.  However,  only 
a  highly  speculative  reason  has  been  put  forward  as  to  why  the  “Heisenberg”  nested  inter¬ 
action  picture  fails  for  most  potentials.  Zhang’s  finite-difference  approach  is  also  an  inter¬ 
esting  option,  requiring  fewer  FFTs  than  the  Lanczos-based  propagation  schemes  tested 
here(46,82).  The  question  of  how  beneficially  the  techniques  developed  here  in  one  dimen¬ 
sion  scale  to  the  much  larger  dimensionalities  of  “real”  molecular  scattering  problems  also 
remains  to  be  settled. 

It  is  important  to  note  that  since  the  short  iterative  Lanczos  propagation  scheme  is 
employed  when  using  the  nested  interaction  picture,  a  greater  number  of  FFTs  is  required 
per  time  step  when  compared  to  the  split-operator  approach  commonly  used  for  time  in¬ 
dependent  Hamiltonians  in  the  SchrQdinger  picture.  This  results  in  a  trade-off  between 
the  computational  savings  afforded  by  the  grid  reduction,  and  the  requirement  for  a  larger 
number  of  FFTs  per  time  step.  An  alternative  to  the  short  iterative  Lanczos  propaga¬ 
tor  that  avoids  this  trade-off  is  the  second  order  finite  difference  propagation  technique 
developed  by  Zhang(46,82).  Using  this  approach,  the  number  of  FFTs  per  time  step  is 
reduced  to  the  same  number  required  by  the  split-operator  method.  This  second  order  fi¬ 
nite  difference  technique  has  been  successfully  applied  to  a  two-dimensional  model  of  CH3I 
photodissociation(83),  and  to  a  three-dimensional  model  of  vibrational  predissociation  of 
van  der  Waals  molecules(84). 

Related  scattering  calculations  into  which  further  investigations  might  be  made  using 
the  interaction  picture  include  time-dependent  potential-energy  surfaces,  molecule-surface 
scattering,  inclusion  of  electronic  degrees  of  freedom  (dropping  the  Born-Oppenheimer.  ap- 
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proximation),  and  adapting  the  models  to  parallel- architecture  computers.  Small,  stochas¬ 
tic,  time-dependent  kicks  to  the  potential  may  be  used  to  model  the  behavior  of  the 
reaction  in  the  presence  of  a  solvent.  The  interaction  picture  is  especially  attractive  for 
dealing  with  explicitly  time-dependent  Hamiltonians,  since  the  Schrodinger-picture  propa¬ 
gator  then  faces  the  same  time-step  issues  that  disadvantage  the  interaction  picture  when 
the  potential  is  time-independent.  The  interaction  of  molecules  with  surfaces  is  of  great 
interest  because  of  the  many  practical  applications  of  such  reactions.  Including  electronic 
degrees  of  freedom  would  make  the  model  more  rigorously  correct,  at  the  cost  of  added 
complexity.  However,  the  Born-Oppenheimer  approximation  is  not  valid  for  all  molecu¬ 
lar  interactions(85).  If  the  prediction  of  equation  (7.3)  holds — in  itself  a  very  important 
question  to  investigate — the  interaction  picture  might  be  employed  to  reduce  the  added 
computational  effort  associated  with  the  addition  of  degrees  of  freedom  associated  with 
atomic  electrons.  Parallelization  is  a  natural  improvement  to  consider  because  the  FFT, 
which  consumes  much  of  the  time  needed  by  these  computations,  has  been  adapted  with 
great  success  to  parallel  computers. 
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Appendix  A.  A  General  Scaling  Rule  for  Computational  Effort  of  Multidimensional 

FFTs 

For  a  one-dimensional  FFT,  the  execution  time  is  governed  by  the  total  number  of  multi¬ 
plications,  which  scales  with  the  grid  size  N  as 

Ei^Nlog^N,  (A.l) 

where  7  is  the  order,  usually  2,  of  the  FFT’s  divide-and-conquer  scheme(86).  Since  we  are 
only  interested  in  relative  scaling  using  the  same  FFT  base,  the  value  of  7  is  unimportant 
and  will  be  dropped  from  now  on. 

Two-dimensional  FFTs  oniVi  x  N2  grids  are  calculated  by  performing  one-dimensional 
FFTs  on  the  Ni  rows  of  length  N2,  followed  by  another  N2  on  the  columns  of  length  Ni. 
Hence,  the  total  effort  scales  with  grid  size  as 

E2  —  N1N2 108  A2  -l-  N2N1  log  Ni 
=  N1N2  log  Ni  N2 

=  P2logP2,  (A.2) 

where  we  define  the  n-dimensional  grid  product, 

Pn  =  flNi.  (A.3) 

i=i 

A  three-dimensional  FFT  breaks  down  into  Ni  two-dimensional  FFTs  on  N2  x  N3 
grids,  plus  N2  on  Ni  x  N3  grids,  and  N3  on  Ni  x  N2  grids.  The  scaling  is  therefore, 

E3  —  Ni  {N2N3  log  N2N3}  +  N2  {N1N3  log  N1N3}  +  N3  {N1N2  log  ^1^2} 

=  PslogCFs)" 

=  2P3logP3.  (A-4) 
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Similarly,  in  four  dimensions  there  are  four  sets  of  three-dimensional  FFTs  involved,  and 
the  scaling  with  grid  size  is. 


Ei  =  Ni  {2JV2iV3iV4  log  N2N3Ni}  +  N2  {2NiNzNi  log  iViAr3iV4} 

+N3  {2iVi  jV2iV4  log  N1N2N4}  +  N4  {2N1N2N3  log  N1N2N3}  (A.5) 

=  2P4log(P4)^ 

=  3!P4logP4.  (A.6) 

The  general  formula  for  computational  effort  of  an  n-dimensional  FFT, 

P„  =  (n-l)!P„logP„,  (A.7) 

is  seen  to  hold  for  n  between  1  and  4.  If  n  >  4,  assuming  (A.7)  to  be  true  for  n  -  1,  the 
scaling  rule  is 

En  =  X^(n-2)!P„log^ 

i=l  " 

=  (n-2)\Pnlogfl^ 

i=i  * 

=  (n-2)!Pnlognj^ 

=  (n-2)\Pnlog^ 

=  (n-2)!PnlogPr' 

=  (n-l)!PnlogP„.  (A.8) 

Equation  (A.7)  is  therefore,  by  induction,  the  scaling  rule  for  all  natural  numbers  n 
of  grid  dimensions. 
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